1
$\begingroup$

I want to fit this function

$$ \operatorname{PF}\left(\sigma_{f_{0}}\right)=1-\exp \left[-\left(\frac{\sigma_{f_{0}}-\sigma_{\min _{0}}}{\lambda}\right)^{k}\right] $$

I already have some points

pt = {
   {73.0, 0.17}, {75.0, 0.33}, {77.0, 0.5}, {79.0, 0.67}, {82.0, 
    0.83}, {83.0, 1}
   };

so I use FindFit function to do this,this is my code

Clear["Global`*"]
(* 数据点 *)
pt = {
   {73.0, 0.17}, {75.0, 0.33}, {77.0, 0.5}, {79.0, 0.67}, {82.0, 
    0.83}, {83.0, 1}
   };

(*点图*)
listPltpt = ListPlot[pt]

(* 基函数 *)
pf[x_] := 1 - Exp[-((x - deltaMin0)/(lambda))^k];
(* 求拟合 *)
ptFit = FindFit[pt, pf@x, {deltaMin0, lambda, k}, x]

(*拟合函数绘图 *)
Show[listPltpt,
 Plot[pf@x /. ptFit, {x, 60, 100}]]

I gave the correct coefficient from the article, and I did the simulation, which is perfect, why can't I figure it out

funC2[del_] := 1 - Exp[-((del - 3.1)/(76.1))^19.8];
pltfunC2 = Plot[funC2@x, {x, 60, 100}];
Show[listPltpt,
 pltfunC2
 ]

enter image description here

=================upgrage========================

=================upgrade========================

by NonlinearModelFit also could not work

enter image description here

The problem was inspired by this excerpt from a paper:

excerpt from paper describing the formula

Since we know that it is a Weibull CDF, then we can just perform a distribution fit:

EstimatedDistribution[
  {73.`, 75.`, 77.`, 79.`, 82.`, 83.`},
  WeibullDistribution[α, β, 3.1],
  WeibullDistribution[19., 76., 3.1]
]
CDF[%][x]

result of CDF

code and results from EstimatedDistribution

Improve this question
$\endgroup$
9
  • $\begingroup$ Have you explored: NonlinearModelFit ? $\endgroup$
    – Syed
    Jun 21 at 9:53
  • $\begingroup$ @Syed I have try,it coud not work $\endgroup$
    – 我心永恒
    Jun 21 at 9:54
  • 1
    $\begingroup$ (1) Please include code that you tried for NonlinearModelFit (2) Update the plot to include the results. Did it improve or did it get worse? (3) Take a look at this post. $\endgroup$
    – Syed
    Jun 21 at 10:11
  • $\begingroup$ @Syed I have upgrage my code and I have see the post $\endgroup$
    – 我心永恒
    Jun 21 at 10:40
  • 3
    $\begingroup$ Would you give a more detailed description of your data? If you have random samples from a Weibull distribution, then you need the counts in each of the bins you seem to have (although the raw observations prior to binning would be even better). In other words the cumulative frequencies of the counts and the boundaries of the bins are essential. Just using the relative frequencies and least squares cannot produce appropriate estimates of precision for the parameters. $\endgroup$
    – JimB
    Jun 21 at 14:39

3 Answers 3

Reset to default
3
$\begingroup$

One way is directly using NMinimize.

Clear["Global`*"]
pt = {{73.0, 0.17}, {75.0, 0.33}, {77.0, 0.5}, {79.0, 0.67}, {82.0, 
    0.83}, {83.0, 1}};
listPltpt = ListPlot[pt];
pf[x_] = 1 - Exp[-((x - deltaMin0)/(lambda))^k];
sol = NMinimize[{Abs[(pf /@ (pt[[;; , 1]]) - pt[[;; , 2]])]^2 // 
    Total, k >= .01}, {deltaMin0, lambda, k},
   Method -> "DifferentialEvolution"]
Show[listPltpt, Plot[pf[x] /. sol[[2]], {x, 70, 90}]]

{0.0084052, {deltaMin0 -> 58.88, lambda -> 19.6469, k -> 4.841}}

enter image description here

Improve this answer
$\endgroup$
3
  • $\begingroup$ It feels good to just do the least square method $\endgroup$
    – 我心永恒
    Jun 21 at 12:11
  • $\begingroup$ Go back to the original method $\endgroup$
    – 我心永恒
    Jun 21 at 12:11
  • $\begingroup$ Only need to set k >= .001. That is k must great than a fixed positive number. $\endgroup$
    – cvgmt
    Jun 21 at 14:36
5
$\begingroup$ 
Clear["Global`*"]

Using NMinimize as shown by cvgmt, you can get good initial estimates for NonlinearModelFit which then provides statistics such as "ANOVATable" or "ParameterTable"

pt = {{73.0, 0.17}, {75.0, 0.33}, {77.0, 0.5}, {79.0, 0.67}, {82.0, 
     0.83}, {83.0, 1}} // Rationalize;

pf[x_] := 1 - Exp[-((x - deltaMin0)/(lambda))^k];

(ptFit = NonlinearModelFit[pt, pf@x, {{deltaMin0, 60}, {lambda, 20}, {k, 5}}, 
      x, WorkingPrecision -> 15])["BestFitParameters"] // N // Quiet

(* {deltaMin0 -> 58.8802, lambda -> 19.6467, k -> 4.84095} *)

ptFit["ANOVATable"] /. x_Real :> N[x]

enter image description here

ptFit["ParameterTable"] /. x_Real :> N[x]

enter image description here

Improve this answer
$\endgroup$
3
  • $\begingroup$ the “BestFitParameters” method ，how can I get this method detail $\endgroup$
    – 我心永恒
    Jun 21 at 15:38
  • $\begingroup$ See the documentation: "NonlinearModelFit returns a symbolic FittedModel object to represent the nonlinear model it constructs. The properties and diagnostics of the model can be obtained from model["property"]." The property "BestFitParameters" gives the parameter estimates from the model. $\endgroup$
    – Bob Hanlon
    Jun 21 at 15:57
  • 1
    $\begingroup$ +1 for including the "ParameterTable". This shows that the data or model or both are inadequate. In addition, ptFit["CorrelationMatrix"] shows near perfect correlations (either +1 or -1) among all 3 parameters. In short, one should not feel very confident about any of the parameter estimates. $\endgroup$
    – JimB
    Jun 22 at 1:46
3
$\begingroup$

It's not just about the numbers but how the numbers were collected. It appears that the OP has digitized Figure 4b in the article mentioned. (I get different values when I digitize that figure but I'll use the numbers obtained by the OP.)

I'm not aware as to how one can directly measure a probability but usually probabilities are estimated with numbers of failures (in this case) out of a fixed number of trials. I'm going to assume that all of the proportions given are the proportion of failures out of 30 trials. The data would then look like

pt = {{73, 30, 5}, {75, 30, 10}, {77, 30, 15}, {79, 30, 19}, {82, 30, 25}, {83, 30, 30}}

So I'm assuming there are 30 trials at each of the $\sigma_{f_0}$ values 73, 75, 77, 79, 82, and 83. The log of the likelihood can be determined as follows.

The probability of failure for a single stress test with $x_i=\sigma_{f_0}(i)$ is

$$p_i=1-\exp \left(-\left(\frac{x_i-\sigma_{Min_0}}{\lambda}\right)^k\right)$$

The log of the likelihood when there are $n_i$ stress tests at $x_i$ with $i=1,2,3,4,5,6$ is

$$\log L=\sum_{i=1}^6 \log\left({{n_i}\choose{y_i}} p_i^{y_i} (1-p_i)^{n_i-y_i}\right)$$

The Mathematica commands are then

pf[x_, deltaMin0_, lambda_, k_] := 1 - Exp[-((x - deltaMin0)/(lambda))^k]

logL = (Log[Binomial[#[[2]], #[[3]]] pf[#[[1]], deltaMin0, lambda, k]^#[[3]]*
  (1 - pf[#[[1]], deltaMin0, lambda, k])^(#[[2]] - #[[3]])] & /@ pt)//Total;

Now the log of the likelihood can be somewhat simplified:

logL = logL //. Log[a_ b_] -> Log[a] + Log[b] //. 
  Log[Exp[a_]] -> a //. Log[a_^b_] -> b Log[a]

The maximum likelihood estimates and associated measures of precision are found as follows:

mle = FindMaximum[{logL, 0 < deltaMin0 < Min[pt[[All, 1]]] && k > 0}, 
  {{deltaMin0, 60}, {lambda, 20}, {k, 5}}]
(* {-11.4408, {deltaMin0 -> 0.00019672, lambda -> 78.7427, k -> 20.8112}} *)

(* Parameter covariance matrix *)
cov = -Inverse[D[logL, {{deltaMin0, lambda, k}, 2}]] /. mle[[2]];
(* Standard errors of estimates *)
se = Sqrt[Diagonal[cov]]
(* {270.198, 270.42, 71.9058} *)
(* Parameter correlation matrix *)
cor = (Table[cov[[i, j]]/Sqrt[cov[[i, i]] cov[[j, j]]], {i, 3}, {j, 3}]) // 
  MatrixForm

Correlation matrix of parameter estimates

The standard errors are huge and the correlations among the parameter estimators are essentially either +1 or -1. If the data is created as I described above, then the results are at best pretty useless if reasonable estimates of all of the parameters are desired. (The resulting estimates are very similar to what is found in the article.)

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.