change point detection of a cummulative distribution

I have a cummulative rainfall time series and I would like to detect the change points. Here's the data.

structure(list(DAY = 1:365, CUMSUM = c(0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3.8, 6.9, 6.9, 6.9, 
6.9, 6.9, 6.9, 6.9, 6.9, 6.9, 6.9, 6.9, 6.9, 6.9, 6.9, 6.9, 7.4, 
7.4, 7.4, 7.4, 7.4, 7.4, 7.4, 7.4, 7.4, 7.4, 7.4, 7.4, 7.4, 22.6, 
22.6, 22.6, 22.6, 22.6, 22.6, 22.8, 26.7, 41.3, 41.3, 44.7, 44.7, 
44.7, 86.8, 92.6, 92.6, 115.2, 117, 126, 134.9, 134.9, 134.9, 
140.7, 140.7, 140.7, 146.5, 146.7, 146.7, 151.7, 152.7, 196.5, 
242.7, 293.4, 331.4, 340, 345.6, 369.5, 442.6, 459, 464.6, 464.6, 
468.2, 475.6, 484.2, 487.8, 487.8, 511, 515, 515, 515, 528.8, 
547.6, 549.4, 549.8, 550, 552.4, 585.9, 798.5, 1062.5, 1107.9, 
1124.5, 1154, 1169.4, 1416.4, 1457.6, 1457.6, 1457.6, 1461.2, 
1464, 1524.7, 1539.5, 1552, 1592.8, 1599.4, 1608.6, 1611.6, 1616.2, 
1656.6, 1667.6, 1667.6, 1668.8, 1680, 1687.1, 1697.9, 1704.7, 
1726.6, 1726.6, 1727.6, 1732.6, 1750.2, 1834.4, 1882.2, 1915.6, 
1940, 1976.6, 2001.2, 2026.4, 2042.6, 2078.1, 2101.2, 2109.2, 
2109.2, 2109.2, 2109.2, 2117, 2117, 2120.2, 2142.4, 2153.4, 2173.4, 
2174.4, 2174.4, 2174.4, 2178.4, 2213.5, 2365.1, 2449.7, 2565.5, 
2673.7, 2749.9, 2830.3, 2896.2, 2920.8, 3236.4, 3266.8, 3288.9, 
3371.5, 3428.5, 3642.5, 3764.9, 3774.9, 3818.7, 3818.7, 3830.9, 
3953.7, 4127.8, 4206, 4217.7, 4217.7, 4219.9, 4220.9, 4220.9, 
4361.1, 4378, 4378, 4388.4, 4393.4, 4417.3, 4419.9, 4419.9, 4419.9, 
4470.3, 4480.3, 4480.7, 4490.7, 4492.9, 4493.4, 4504, 4504, 4504, 
4505.4, 4509.8, 4509.8, 4509.8, 4509.8, 4509.8, 4509.8, 4509.8, 
4510.4, 4510.4, 4512.8, 4515.4, 4517.8, 4527.5, 4532.1, 4539.7, 
4541.7, 4573.3, 4606.5, 4607.3, 4613.5, 4613.5, 4613.5, 4613.5, 
4613.5, 4613.5, 4613.5, 4613.5, 4613.5, 4613.5, 4613.9, 4621.1, 
4621.1, 4621.1, 4636.5, 4647.9, 4649.1, 4649.3, 4649.3, 4649.3, 
4655, 4655, 4663.6, 4663.6, 4664.2, 4664.2, 4665, 4665, 4665, 
4665, 4665, 4665, 4665, 4665, 4665, 4665, 4665, 4665, 4665, 4665, 
4665.9, 4665.9, 4665.9, 4665.9, 4665.9, 4665.9, 4665.9, 4665.9, 
4665.9, 4665.9, 4665.9, 4665.9, 4665.9, 4673.1, 4673.1, 4673.1, 
4673.1, 4673.1, 4673.1, 4673.1, 4673.1, 4673.1, 4673.5, 4673.5, 
4673.5, 4673.5, 4673.5, 4673.5, 4673.5, 4673.5, 4673.5)), .Names = 
c("DAY","CUMSUM"), class = "data.frame", row.names = c(NA, -365L))

I would like to apply a two phase linear regression in detecting the change points here using R.

There is a matlab code available here https://www.mathworks.com/matlabcentral/fileexchange/26804-two-phase-linear-regression-model

but there is no equivalent package in R.

Can anyone suggest how to do this?

Here's the expected output.

We can use the R package segmented; here is a step-by-step example.

1. Load the library.

   library(segmented);
   

2. Fit a piecewise linear model with two breakpoints to the sample data (here I assume df contains the data as a data.frame). Note that we must provide some guesses for the breakpoints.

   fit <- lm(CUMSUM ~ DAY, data = df);
   fit.seg <- segmented(fit, psi = c(100, 200));
   fit.seg;
   #Call: segmented.lm(obj = fit, psi = c(100, 200))
   #
   #Meaningful coefficients of the linear terms:
   #(Intercept)          DAY       U1.DAY       U2.DAY
   #     -58.20         1.25        35.70       -34.98
   #
   #Estimated Break-Point(s):
   #psi1.DAY  psi2.DAY
   #   153.8     272.9
   

3. We plot the curve and mark the breakpoint estimates in red.

   library(ggplot2);
   ggplot(df, aes(DAY, CUMSUM)) +
       geom_line() +
       geom_vline(data = as.data.frame(fit.seg$psi), aes(xintercept = `Est.`), col = "red")
   

4. More details can be found in the segmented reference manual on CRAN. The return object fit.seg also contains parameter estimates for each piece.