solving a equation with integral

1 vue (au cours des 30 derniers jours)
Louis Liu
Louis Liu le 22 Août 2017
Modifié(e) : Louis Liu le 29 Août 2017
Hello,
I write some code like this,
>> C = 1.00;
>> n = 38;
>> alpha = 0.05;
>> Cp = C+0.33;
>> fun=@(y)chi2cdf((n-1)*(3*Cp*sqrt(n)-y).^2/(9*n*czero^2),n-1).*(normpdf(y+3*(Cp-C)*sqrt(n))+normpdf(y-3*(Cp-C)*sqrt(n)));
>> solve(integral(fun,0,3*Cp*sqrt(n))-alpha,czero)
and get the main error message : Undefined function or variable 'czero'.
Does anyone can tell me how to solve this equation? Using symbolic variable or other method?
Thanks!

Réponse acceptée

Walter Roberson
Walter Roberson le 22 Août 2017
What you would like to do is roughly
C = 1.00;
n = 38;
alpha = 0.05;
Cp = C+0.33;
syms y czero real
fun = chi2cdf((n-1)*(3*Cp*sqrt(n)-y).^2/(9*n*czero^2),n-1).*(normpdf(y+3*(Cp-C)*sqrt(n))+normpdf(y-3*(Cp-C)*sqrt(n)));
solve( int(fun, y, 0,3*Cp*sqrt(n))-alpha, czero)
Unfortunately, chi2cdf does not expect symbolic values, and naively does a test on whether the input < 0 (because input values below 0 give an output of 0), but that test does not work on symbolic values. You can avoid the test by coding chi2pdf in terms of gampdf() but that only postpones the problem, as gampdf() has the same difficulty.
So, instead you have to work numerically:
C = 1.00;
n = 38;
alpha = 0.05;
Cp = C+0.33;
fun = @(y, czero) chi2cdf((n-1)*(3*Cp*sqrt(n)-y).^2/(9*n*czero^2),n-1).*(normpdf(y+3*(Cp-C)*sqrt(n))+normpdf(y-3*(Cp-C)*sqrt(n)));
czero_guess = 0;
czero_sol = fzero( @(czero) integral( @(y) fun(y, czero), 0, 3*Cp*sqrt(n) ) - alpha, czero_guess )
The output I am getting is -1.26021065783808
  3 commentaires
Walter Roberson
Walter Roberson le 22 Août 2017
I had to choose something for the guess, and I had no idea of the value that would be found.
The only place your equation has czero is in the form czero^2 so your equation cannot tell the difference between negative and positive czero
Louis Liu
Louis Liu le 23 Août 2017
Oh, you are right! Thanks!

Connectez-vous pour commenter.

Plus de réponses (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!