could someone tell me how to implement a multiple comparison test for proportions in a 2xc crosstabulation

classic Classic list List threaded Threaded
17 messages Options
Reply | Threaded
Open this post in threaded view
|

could someone tell me how to implement a multiple comparison test for proportions in a 2xc crosstabulation

xingwang ye
Dear all,
I wanna to do multiple comparison test for proportions (multiple chi
squre ?), could someone tell me how in R, thank you!

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Re: could someone tell me how to implement a multiple comparison test for proportions in a 2xc crosstabulation

Wincent
?p.adjust
It's the function used to Adjust P-values for Multiple Comparisons.


2006/6/14, xingwang ye <[hidden email]>:
> Dear all,
> I wanna to do multiple comparison test for proportions (multiple chi
> squre ?), could someone tell me how in R, thank you!
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>


--
黄荣贵
Deparment of Sociology
Fudan University


______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Re: could someone tell me how to implement a multiple comparison test for proportions in a 2xc crosstabulation

Xiaoting Hua
In reply to this post by xingwang ye
Maybe you want the package "multcomp".

On 6/14/06, xingwang ye <[hidden email]> wrote:
> Dear all,
> I wanna to do multiple comparison test for proportions (multiple chi
> squre ?), could someone tell me how in R, thank you!
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>


--
  此致
敬礼!
  华孝挺
浙江大学核农所
杭州,中国
310029


______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

number of iteration s exceeded maximum of 50

Leaf Sun
In reply to this post by Wincent
Hi all,

I found r-site-research not work for me these days.

When I was doing nls( ) , there was an error "number of iterations exceeded maximum of 50". I set number in nls.control which is supposed to control the number of iterations but it didn't work well. Could anybody with this experience tell me how to fix it? Thanks in advance!

Leaf

        [[alternative HTML version deleted]]

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Re: number of iteration s exceeded maximum of 50

Uwe Ligges
Leaf Sun wrote:

> Hi all,
>
> I found r-site-research not work for me these days.
>
> When I was doing nls( ) , there was an error "number of iterations exceeded maximum of 50". I set number in nls.control which is supposed to control the number of iterations but it didn't work well. Could anybody with this experience tell me how to fix it? Thanks in advance!

We cannot make suggestions unless you tell us what you tried yourself.
Id possible, please gib´ve a reproducible examle.

Uwe Ligges

> Leaf
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Re: number of iteration s exceeded maximum of 50

Leaf Sun
Sorry, I thought it was a straightforward question inside which I was stuck .

I used nls( ) to estimate a and b in this function.

nls(y~ a*x^b,start=list(a=a1,b=b1)

seems the start list I gave was not able to reach convergence and it gave notes: number of iteration s exceeded maximum of 50. Then I put  nls.control(maxiter = 50, tol = 1e-05, minFactor = 1/1024) in nls(.. ), and modified the argument of maxiter = 500. But it worked out as the same way and noted : number of iteration s exceeded maximum of 50. I have totally no idea how to set this parameter MAXITER.

Thanks for any information!

Leaf


>  Hi  all,
>  
>  I  found  r-site-research  not  work  for  me  these  days.
>  
>  When  I  was  doing  nls(  )  ,  there  was  an  error  "number  of  iterations  exceeded  maximum  of  50".  I  set  number  in  nls.control  which  is  supposed  to  control  the  number  of  iterations  but  it  didn't  work  well.  Could  anybody  with  this  experience  tell  me  how  to  fix  it?  Thanks  in  advance!

We  cannot  make  suggestions  unless  you  tell  us  what  you  tried  yourself.
Id  possible,  please  gib´ve  a  reproducible  examle.

Uwe  Ligges

>  Leaf
>  
>   [[alternative  HTML  version  deleted]]
>  
>  ______________________________________________
>  [hidden email]  mailing  list
>  https://stat.ethz.ch/mailman/listinfo/r-help
>  PLEASE  do  read  the  posting  guide!  http://www.R-project.org/posting-guide.html

        [[alternative HTML version deleted]]


______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Re: number of iteration s exceeded maximum of 50

Bert Gunter
My goodness! This is **NOT** a reproducible example. You need to give us the
exact data you fitted to reproduce your results/diagnose your problem.


-- Bert Gunter
Genentech Non-Clinical Statistics
South San Francisco, CA
 
"The business of the statistician is to catalyze the scientific learning
process."  - George E. P. Box
 
 

> -----Original Message-----
> From: [hidden email]
> [mailto:[hidden email]] On Behalf Of Leaf Sun
> Sent: Friday, June 16, 2006 9:41 AM
> To: Uwe Ligges
> Cc: [hidden email]
> Subject: Re: [R] number of iteration s exceeded maximum of 50
>
> Sorry, I thought it was a straightforward question inside
> which I was stuck .
>
> I used nls( ) to estimate a and b in this function.
>
> nls(y~ a*x^b,start=list(a=a1,b=b1)
>
> seems the start list I gave was not able to reach convergence
> and it gave notes: number of iteration s exceeded maximum of
> 50. Then I put  nls.control(maxiter = 50, tol = 1e-05,
> minFactor = 1/1024) in nls(.. ), and modified the argument of
> maxiter = 500. But it worked out as the same way and noted :
> number of iteration s exceeded maximum of 50. I have totally
> no idea how to set this parameter MAXITER.
>
> Thanks for any information!
>
> Leaf
>
>
> >  Hi  all,
> >  
> >  I  found  r-site-research  not  work  for  me  these  days.
> >  
> >  When  I  was  doing  nls(  )  ,  there  was  an  error  
> "number  of  iterations  exceeded  maximum  of  50".  I  set  
> number  in  nls.control  which  is  supposed  to  control  
> the  number  of  iterations  but  it  didn't  work  well.  
> Could  anybody  with  this  experience  tell  me  how  to  
> fix  it?  Thanks  in  advance!
>
> We  cannot  make  suggestions  unless  you  tell  us  what  
> you  tried  yourself.
> Id  possible,  please  gib4ve  a  reproducible  examle.
>
> Uwe  Ligges
>
> >  Leaf
> >  
> >   [[alternative  HTML  version  deleted]]
> >  
> >  ______________________________________________
> >  [hidden email]  mailing  list
> >  https://stat.ethz.ch/mailman/listinfo/r-help
> >  PLEASE  do  read  the  posting  guide!  
> http://www.R-project.org/posting-guide.html
>
> [[alternative HTML version deleted]]
>
>

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Re: number of iteration s exceeded maximum of 50

Douglas Bates-2
In reply to this post by Leaf Sun
On 6/16/06, Leaf Sun <[hidden email]> wrote:
> Sorry, I thought it was a straightforward question inside which I was stuck .
>
> I used nls( ) to estimate a and b in this function.
>
> nls(y~ a*x^b,start=list(a=a1,b=b1)
>
> seems the start list I gave was not able to reach convergence and it gave notes: number of iteration s exceeded maximum of 50. Then I put  nls.control(maxiter = 50, tol = 1e-05, minFactor = 1/1024) in nls(.. ), and modified the argument of maxiter = 500. But it worked out as the same way and noted : number of iteration s exceeded maximum of 50. I have totally no idea how to set this parameter MAXITER.
>
> Thanks for any information!

I think you are assuming that values passed to nls.control are
persistent and will apply to further calls to nls.  They don't.  If
you want to increase the maximum number of iterations you do it as

nls(y ~ a*x^b, start = list(a = a1, b = b1), control = list(maxiter = 500))

but I would suggest that you also use trace = TRUE in the call to nls
so you can see where the iterations are going.  Merely increasing the
number of iterations for an optimization that has gone into
never-never land isn't going to help it converge.

Two other things to consider: this is a partially linear model in the
the parameter `a' appears linearly in the model expression.  You may
be able to stabilize the iterations using

 nls(y ~ x^b, start = list(b = b1), control = list(maxiter = 500),
trace = TRUE, alg = 'plinear')

Finally, and most important, please plot the data before trying to fit
a nonlinear model to it so you can see if it has the characteristics
that you would expect from data generate by such a model.  As Brian
Joiner said, "Regression without plots is truly a regression".

>
> Leaf
>
>
> >  Hi  all,
> >
> >  I  found  r-site-research  not  work  for  me  these  days.
> >
> >  When  I  was  doing  nls(  )  ,  there  was  an  error  "number  of  iterations  exceeded  maximum  of  50".  I  set  number  in  nls.control  which  is  supposed  to  control  the  number  of  iterations  but  it  didn't  work  well.  Could  anybody  with  this  experience  tell  me  how  to  fix  it?  Thanks  in  advance!
>
> We  cannot  make  suggestions  unless  you  tell  us  what  you  tried  yourself.
> Id  possible,  please  gib´ve  a  reproducible  examle.
>
> Uwe  Ligges
>
> >  Leaf
> >
> >   [[alternative  HTML  version  deleted]]
> >
> >  ______________________________________________
> >  [hidden email]  mailing  list
> >  https://stat.ethz.ch/mailman/listinfo/r-help
> >  PLEASE  do  read  the  posting  guide!  http://www.R-project.org/posting-guide.html
>
>         [[alternative HTML version deleted]]
>
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>
>

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Re: number of iteration s exceeded maximum of 50

Leaf Sun
Thanks to Douglas and all others who responded.

I applied  nls(y  ~  a*x^b,  start  =  list(a  =  a1,  b  =  b1),  control  =  list(maxiter  =  500), trace=TRUE) to increase the number of iterations, found it successful. The suggestion Douglas raised in plotting the data and then tracing the optim numbers is correct because I found when I gave the number of b1 oppositely(say, should be positive, then given negative), nls( ) would never reached the convergence. Thanks for the nice suggestions!

Leaf

>  Sorry,  I  thought  it  was  a  straightforward  question  inside  which  I  was  stuck  .
>
>  I  used  nls(  )  to  estimate  a  and  b  in  this  function.
>
>  nls(y~  a*x^b,start=list(a=a1,b=b1)
>
>  seems  the  start  list  I  gave  was  not  able  to  reach  convergence  and  it  gave  notes:  number  of  iteration  s  exceeded  maximum  of  50.  Then  I  put    nls.control(maxiter  =  50,  tol  =  1e-05,  minFactor  =  1/1024)  in  nls(..  ),  and  modified  the  argument  of  maxiter  =  500.  But  it  worked  out  as  the  same  way  and  noted  :  number  of  iteration  s  exceeded  maximum  of  50.  I  have  totally  no  idea  how  to  set  this  parameter  MAXITER.
>
>  Thanks  for  any  information!

I  think  you  are  assuming  that  values  passed  to  nls.control  are
persistent  and  will  apply  to  further  calls  to  nls.    They  don't.    If
you  want  to  increase  the  maximum  number  of  iterations  you  do  it  as

nls(y  ~  a*x^b,  start  =  list(a  =  a1,  b  =  b1),  control  =  list(maxiter  =  500))

but  I  would  suggest  that  you  also  use  trace  =  TRUE  in  the  call  to  nls
so  you  can  see  where  the  iterations  are  going.    Merely  increasing  the
number  of  iterations  for  an  optimization  that  has  gone  into
never-never  land  isn't  going  to  help  it  converge.

Two  other  things  to  consider:  this  is  a  partially  linear  model  in  the
the  parameter  `a'  appears  linearly  in  the  model  expression.    You  may
be  able  to  stabilize  the  iterations  using

 nls(y  ~  x^b,  start  =  list(b  =  b1),  control  =  list(maxiter  =  500),
trace  =  TRUE,  alg  =  'plinear')

Finally,  and  most  important,  please  plot  the  data  before  trying  to  fit
a  nonlinear  model  to  it  so  you  can  see  if  it  has  the  characteristics
that  you  would  expect  from  data  generate  by  such  a  model.    As  Brian
Joiner  said,  "Regression  without  plots  is  truly  a  regression".

>
>  Leaf
>
>
>   >    Hi    all,
>   >
>   >    I    found    r-site-research    not    work    for    me    these    days.
>   >
>   >    When    I    was    doing    nls(    )    ,    there    was    an    error    "number    of    iterations    exceeded    maximum    of    50".    I    set    number    in    nls.control    which    is    supposed    to    control    the    number    of    iterations    but    it    didn't    work    well.    Could    anybody    with    this    experience    tell    me    how    to    fix    it?    Thanks    in    advance!
>
>  We    cannot    make    suggestions    unless    you    tell    us    what    you    tried    yourself.
>  Id    possible,    please    gib´ve    a    reproducible    examle.
>
>  Uwe    Ligges
>
>   >    Leaf
>   >
>   >      [[alternative    HTML    version    deleted]]
>   >
>   >    ______________________________________________
>   >    [hidden email]    mailing    list
>   >    https://stat.ethz.ch/mailman/listinfo/r-help
>   >    PLEASE    do    read    the    posting    guide!    http://www.R-project.org/posting-guide.html
>
>                  [[alternative  HTML  version  deleted]]
>
>
>  ______________________________________________
>  [hidden email]  mailing  list
>  https://stat.ethz.ch/mailman/listinfo/r-help
>  PLEASE  do  read  the  posting  guide!  http://www.R-project.org/posting-guide.html
>
>
        [[alternative HTML version deleted]]


______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Weibull distribution

Leaf Sun
Hi all,

By its definition, the mean and variance of two-par. Weibull distribution are:

 

 

 (www.wikipedia.org)


I was wondering, if given mean and sd. could we parameterize the distribution? I tried this in R.

gamma.fun <- function(mu,sd,start=100)    
{
f.fn <- function(alpha) sd^2-mu^2/(gamma(1+1/alpha))^2*(gamma(1+2/alpha)-(gamma(1+1/alpha))^2)
alpha <- optim(start, f.fn,method='BFGS')
beta <- mu/gamma(1+1/alpha$par)
return(list=c(a=alpha$par,b=beta));
}


But the problems come up here:

1)  the return values of a and b are only related to the input mean, and nothing to do with the sd. For instance, when I apply a mean mu = 3 whatever I use sd=2, sd=4, the function returned the same scale and shape values.

> gamma.fun(3,4,10);
       a        b
5.112554 3.263178

> gamma.fun(3,2,10);
       a        b
5.112554 3.263178

2) the start value determines the results: if I apply mean = 3, and sd=2, with a start of 10, it would return alpha close to 10, if I use a start = 100, it would return alpha close to 100.

> gamma.fun(3,2,10);
       a        b
5.112554 3.263178

> gamma.fun(3,2,100);
        a         b
99.999971  3.017120

Since I am not a statistician, I guess there must be some theoretical reasons wrong with this question. So I am looking forward to some correction and advice to solve these. Thanks a lot in advance!

Leaf

        [[alternative HTML version deleted]]

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Re: Weibull distribution

Leaf Sun
Hi William,

Thanks a lot for your response. I checked the package and found that what I want to solve was the opposite, that is, from mean and sd to parameters shape and scale. Could anyone give some hints please? Any suggestion would be appreciated!

Leaf



----- Original Message -----

From: William Asquith,  [hidden email]
Sent: 2006-07-17,  16:18:31
To: Leaf Sun, [hidden email]
Subject:  Re: [R] Weibull distribution
 
Do  not  have  answer  per  se,  but  if  you  are  seeking  some  comparisons--  
try  three  parameter  Weibull  as  implemented  by  the  lmomco  package.

William
On  Jul  17,  2006,  at  1:18  PM,  Leaf  Sun  wrote:

>  Hi  all,
>
>  By  its  definition,  the  mean  and  variance  of  two-par.  Weibull    
>  distribution  are:
>
>
>
>
>
>    (www.wikipedia.org)
>
>
>  I  was  wondering,  if  given  mean  and  sd.  could  we  parameterize  the    
>  distribution?  I  tried  this  in  R.
>
>  gamma.fun   <-  function(mu,sd,start=100)
>  {
>  f.fn   <-  function(alpha)  sd^2-mu^2/(gamma(1+1/alpha))^2*(gamma(1+2/  
>  alpha)-(gamma(1+1/alpha))^2)
>  alpha   <-  optim(start,  f.fn,method='BFGS')
>  beta   <-  mu/gamma(1+1/alpha$par)
>  return(list=c(a=alpha$par,b=beta));
>  }
>
>
>  But  the  problems  come  up  here:
>
>  1)    the  return  values  of  a  and  b  are  only  related  to  the  input    
>  mean,  and  nothing  to  do  with  the  sd.  For  instance,  when  I  apply  a    
>  mean  mu  =  3  whatever  I  use  sd=2,  sd=4,  the  function  returned  the    
>  same  scale  and  shape  values.
>
> >  gamma.fun(3,4,10);
>                a                b
>  5.112554  3.263178
>
> >  gamma.fun(3,2,10);
>                a                b
>  5.112554  3.263178
>
>  2)  the  start  value  determines  the  results:  if  I  apply  mean  =  3,  and    
>  sd=2,  with  a  start  of  10,  it  would  return  alpha  close  to  10,  if  I    
>  use  a  start  =  100,  it  would  return  alpha  close  to  100.
>
> >  gamma.fun(3,2,10);
>                a                b
>  5.112554  3.263178
>
> >  gamma.fun(3,2,100);
>                  a                  b
>  99.999971    3.017120
>
>  Since  I  am  not  a  statistician,  I  guess  there  must  be  some    
>  theoretical  reasons  wrong  with  this  question.  So  I  am  looking    
>  forward  to  some  correction  and  advice  to  solve  these.  Thanks  a  lot    
>  in  advance!
>
>  Leaf
>
>   [[alternative  HTML  version  deleted]]
>
>  ______________________________________________
>  [hidden email]  mailing  list
>  https://stat.ethz.ch/mailman/listinfo/r-help
>  PLEASE  do  read  the  posting  guide!  http://www.R-project.org/posting- 
>  guide.html

        [[alternative HTML version deleted]]

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

(no subject)

Leaf Sun
In reply to this post by Leaf Sun
@yahoo.ca>
Subject:  Re: [R] Weibull distribution
Message-ID: <[hidden email]>
X-mailer: Foxmail 6, 3, 103, 21 [cn]
Mime-Version: 1.0
Content-Type: multipart/alternative;
        boundary="=====003_Dragon527446281311_====="


This is a multi-part message in MIME format.

--=====003_Dragon527446281311_=====
Content-Type: text/plain;
        charset="gb2312"
Content-Transfer-Encoding: 7bit

Hi William,

Thanks a lot for your response. I checked the package and found that what I want to solve was the opposite, that is, from mean and sd to parameters shape and scale. Could anyone give some hints please? Any suggestion would be appreciated!

Leaf



----- Original Message -----

From: William Asquith,  [hidden email]
Sent: 2006-07-17,  16:18:31
To: Leaf Sun, [hidden email]
Subject:  Re: [R] Weibull distribution
 
Do  not  have  answer  per  se,  but  if  you  are  seeking  some  comparisons--  
try  three  parameter  Weibull  as  implemented  by  the  lmomco  package.

William
On  Jul  17,  2006,  at  1:18  PM,  Leaf  Sun  wrote:

>  Hi  all,
>
>  By  its  definition,  the  mean  and  variance  of  two-par.  Weibull    
>  distribution  are:
>
>
>
>
>
>    (www.wikipedia.org)
>
>
>  I  was  wondering,  if  given  mean  and  sd.  could  we  parameterize  the    
>  distribution?  I  tried  this  in  R.
>
>  gamma.fun   <-  function(mu,sd,start=100)
>  {
>  f.fn   <-  function(alpha)  sd^2-mu^2/(gamma(1+1/alpha))^2*(gamma(1+2/  
>  alpha)-(gamma(1+1/alpha))^2)
>  alpha   <-  optim(start,  f.fn,method='BFGS')
>  beta   <-  mu/gamma(1+1/alpha$par)
>  return(list=c(a=alpha$par,b=beta));
>  }
>
>
>  But  the  problems  come  up  here:
>
>  1)    the  return  values  of  a  and  b  are  only  related  to  the  input    
>  mean,  and  nothing  to  do  with  the  sd.  For  instance,  when  I  apply  a    
>  mean  mu  =  3  whatever  I  use  sd=2,  sd=4,  the  function  returned  the    
>  same  scale  and  shape  values.
>
> >  gamma.fun(3,4,10);
>                a                b
>  5.112554  3.263178
>
> >  gamma.fun(3,2,10);
>                a                b
>  5.112554  3.263178
>
>  2)  the  start  value  determines  the  results:  if  I  apply  mean  =  3,  and    
>  sd=2,  with  a  start  of  10,  it  would  return  alpha  close  to  10,  if  I    
>  use  a  start  =  100,  it  would  return  alpha  close  to  100.
>
> >  gamma.fun(3,2,10);
>                a                b
>  5.112554  3.263178
>
> >  gamma.fun(3,2,100);
>                  a                  b
>  99.999971    3.017120
>
>  Since  I  am  not  a  statistician,  I  guess  there  must  be  some    
>  theoretical  reasons  wrong  with  this  question.  So  I  am  looking    
>  forward  to  some  correction  and  advice  to  solve  these.  Thanks  a  lot    
>  in  advance!
>
>  Leaf
>
>   [[alternative  HTML  version  deleted]]
>
>  ______________________________________________
>  [hidden email]  mailing  list
>  https://stat.ethz.ch/mailman/listinfo/r-help
>  PLEASE  do  read  the  posting  guide!  http://www.R-project.org/posting- 
>  guide.html

--=====003_Dragon527446281311_=====
Content-Type: text/html;
        charset="gb2312"
Content-Transfer-Encoding: 7bit

<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<HTML><HEAD>
<META http-equiv=Content-Type content="text/html; charset=gb2312">
<META content="MSHTML 6.00.2900.2912" name=GENERATOR></HEAD>
<BODY>
<DIV>Hi William,</DIV>
<DIV>&nbsp;</DIV>
<DIV>Thanks a lot for your response. I checked the package and found that what I
want to solve was the opposite, that is, from mean and sd to parameters shape
and scale. Could anyone give some hints please? Any suggestion would be
appreciated!</DIV>
<DIV><BR>Leaf</DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;</DIV>
<DIV>----- Original Message -----</DIV>
<DIV>&nbsp;</DIV>
<DIV><FONT size=2><FONT face=Tahoma><STRONG>From:</STRONG> William
Asquith,&nbsp;&nbsp;<A
href="mailto:[hidden email]">[hidden email]</A><BR><B>Sent:</B>
2006-07-17,&nbsp; 16:18:31<BR><B>To:</B> Leaf Sun, <A
href="mailto:[hidden email]">[hidden email]</A><BR><B>Subject:</B>&nbsp;
Re: [R] Weibull distribution</FONT></FONT></DIV>
<DIV>&nbsp;&nbsp;</DIV>
<DIV>
<TABLE width="100%">
  <TBODY>
  <TR>
    <TD width="100%">
      <BLOCKQUOTE
      style="PADDING-RIGHT: 0px; PADDING-LEFT: 5px; MARGIN-LEFT: 5px; BORDER-LEFT: #000000 2px solid; MARGIN-RIGHT: 0px">
        <DIV>Do &nbsp;not &nbsp;have &nbsp;answer &nbsp;per &nbsp;se, &nbsp;but
        &nbsp;if &nbsp;you &nbsp;are &nbsp;seeking &nbsp;some
        &nbsp;comparisons-- &nbsp;</DIV>
        <DIV>try &nbsp;three &nbsp;parameter &nbsp;Weibull &nbsp;as
        &nbsp;implemented &nbsp;by &nbsp;the &nbsp;lmomco &nbsp;package.</DIV>
        <DIV>&nbsp;</DIV>
        <DIV>William</DIV>
        <DIV>On &nbsp;Jul &nbsp;17, &nbsp;2006, &nbsp;at &nbsp;1:18 &nbsp;PM,
        &nbsp;Leaf &nbsp;Sun &nbsp;wrote:</DIV>
        <DIV>&nbsp;</DIV>
        <DIV>&gt; &nbsp;Hi &nbsp;all,</DIV>
        <DIV>&gt;</DIV>
        <DIV>&gt; &nbsp;By &nbsp;its &nbsp;definition, &nbsp;the &nbsp;mean
        &nbsp;and &nbsp;variance &nbsp;of &nbsp;two-par. &nbsp;Weibull &nbsp;
        &nbsp;</DIV>
        <DIV>&gt; &nbsp;distribution &nbsp;are:</DIV>
        <DIV>&gt;</DIV>
        <DIV>&gt;</DIV>
        <DIV>&gt;</DIV>
        <DIV>&gt;</DIV>
        <DIV>&gt;</DIV>
        <DIV>&gt; &nbsp; &nbsp;(www.wikipedia.org)</DIV>
        <DIV>&gt;</DIV>
        <DIV>&gt;</DIV>
        <DIV>&gt; &nbsp;I &nbsp;was &nbsp;wondering, &nbsp;if &nbsp;given
        &nbsp;mean &nbsp;and &nbsp;sd. &nbsp;could &nbsp;we &nbsp;parameterize
        &nbsp;the &nbsp; &nbsp;</DIV>
        <DIV>&gt; &nbsp;distribution? &nbsp;I &nbsp;tried &nbsp;this &nbsp;in
        &nbsp;R.</DIV>
        <DIV>&gt;</DIV>
        <DIV>&gt; &nbsp;gamma.fun &nbsp; &lt;-
        &nbsp;function(mu,sd,start=100)</DIV>
        <DIV>&gt; &nbsp;{</DIV>
        <DIV>&gt; &nbsp;f.fn &nbsp; &lt;- &nbsp;function(alpha)
        &nbsp;sd^2-mu^2/(gamma(1+1/alpha))^2*(gamma(1+2/ &nbsp;</DIV>
        <DIV>&gt; &nbsp;alpha)-(gamma(1+1/alpha))^2)</DIV>
        <DIV>&gt; &nbsp;alpha &nbsp; &lt;- &nbsp;optim(start,
        &nbsp;f.fn,method='BFGS')</DIV>
        <DIV>&gt; &nbsp;beta &nbsp; &lt;- &nbsp;mu/gamma(1+1/alpha$par)</DIV>
        <DIV>&gt; &nbsp;return(list=c(a=alpha$par,b=beta));</DIV>
        <DIV>&gt; &nbsp;}</DIV>
        <DIV>&gt;</DIV>
        <DIV>&gt;</DIV>
        <DIV>&gt; &nbsp;But &nbsp;the &nbsp;problems &nbsp;come &nbsp;up
        &nbsp;here:</DIV>
        <DIV>&gt;</DIV>
        <DIV>&gt; &nbsp;1) &nbsp; &nbsp;the &nbsp;return &nbsp;values &nbsp;of
        &nbsp;a &nbsp;and &nbsp;b &nbsp;are &nbsp;only &nbsp;related &nbsp;to
        &nbsp;the &nbsp;input &nbsp; &nbsp;</DIV>
        <DIV>&gt; &nbsp;mean, &nbsp;and &nbsp;nothing &nbsp;to &nbsp;do
        &nbsp;with &nbsp;the &nbsp;sd. &nbsp;For &nbsp;instance, &nbsp;when
        &nbsp;I &nbsp;apply &nbsp;a &nbsp; &nbsp;</DIV>
        <DIV>&gt; &nbsp;mean &nbsp;mu &nbsp;= &nbsp;3 &nbsp;whatever &nbsp;I
        &nbsp;use &nbsp;sd=2, &nbsp;sd=4, &nbsp;the &nbsp;function
        &nbsp;returned &nbsp;the &nbsp; &nbsp;</DIV>
        <DIV>&gt; &nbsp;same &nbsp;scale &nbsp;and &nbsp;shape
        &nbsp;values.</DIV>
        <DIV>&gt;</DIV>
        <DIV>&gt; &gt; &nbsp;gamma.fun(3,4,10);</DIV>
        <DIV>&gt; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;a
        &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;b</DIV>
        <DIV>&gt; &nbsp;5.112554 &nbsp;3.263178</DIV>
        <DIV>&gt;</DIV>
        <DIV>&gt; &gt; &nbsp;gamma.fun(3,2,10);</DIV>
        <DIV>&gt; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;a
        &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;b</DIV>
        <DIV>&gt; &nbsp;5.112554 &nbsp;3.263178</DIV>
        <DIV>&gt;</DIV>
        <DIV>&gt; &nbsp;2) &nbsp;the &nbsp;start &nbsp;value &nbsp;determines
        &nbsp;the &nbsp;results: &nbsp;if &nbsp;I &nbsp;apply &nbsp;mean &nbsp;=
        &nbsp;3, &nbsp;and &nbsp; &nbsp;</DIV>
        <DIV>&gt; &nbsp;sd=2, &nbsp;with &nbsp;a &nbsp;start &nbsp;of &nbsp;10,
        &nbsp;it &nbsp;would &nbsp;return &nbsp;alpha &nbsp;close &nbsp;to
        &nbsp;10, &nbsp;if &nbsp;I &nbsp; &nbsp;</DIV>
        <DIV>&gt; &nbsp;use &nbsp;a &nbsp;start &nbsp;= &nbsp;100, &nbsp;it
        &nbsp;would &nbsp;return &nbsp;alpha &nbsp;close &nbsp;to
        &nbsp;100.</DIV>
        <DIV>&gt;</DIV>
        <DIV>&gt; &gt; &nbsp;gamma.fun(3,2,10);</DIV>
        <DIV>&gt; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;a
        &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;b</DIV>
        <DIV>&gt; &nbsp;5.112554 &nbsp;3.263178</DIV>
        <DIV>&gt;</DIV>
        <DIV>&gt; &gt; &nbsp;gamma.fun(3,2,100);</DIV>
        <DIV>&gt; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;
        &nbsp;a &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;
        &nbsp;b</DIV>
        <DIV>&gt; &nbsp;99.999971 &nbsp; &nbsp;3.017120</DIV>
        <DIV>&gt;</DIV>
        <DIV>&gt; &nbsp;Since &nbsp;I &nbsp;am &nbsp;not &nbsp;a
        &nbsp;statistician, &nbsp;I &nbsp;guess &nbsp;there &nbsp;must &nbsp;be
        &nbsp;some &nbsp; &nbsp;</DIV>
        <DIV>&gt; &nbsp;theoretical &nbsp;reasons &nbsp;wrong &nbsp;with
        &nbsp;this &nbsp;question. &nbsp;So &nbsp;I &nbsp;am &nbsp;looking
        &nbsp; &nbsp;</DIV>
        <DIV>&gt; &nbsp;forward &nbsp;to &nbsp;some &nbsp;correction &nbsp;and
        &nbsp;advice &nbsp;to &nbsp;solve &nbsp;these. &nbsp;Thanks &nbsp;a
        &nbsp;lot &nbsp; &nbsp;</DIV>
        <DIV>&gt; &nbsp;in &nbsp;advance!</DIV>
        <DIV>&gt;</DIV>
        <DIV>&gt; &nbsp;Leaf</DIV>
        <DIV>&gt;</DIV>
        <DIV>&gt; &nbsp; [[alternative &nbsp;HTML &nbsp;version
        &nbsp;deleted]]</DIV>
        <DIV>&gt;</DIV>
        <DIV>&gt; &nbsp;______________________________________________</DIV>
        <DIV>&gt; &nbsp;[hidden email] &nbsp;mailing &nbsp;list</DIV>
        <DIV>&gt; &nbsp;https://stat.ethz.ch/mailman/listinfo/r-help</DIV>
        <DIV>&gt; &nbsp;PLEASE &nbsp;do &nbsp;read &nbsp;the &nbsp;posting
        &nbsp;guide! &nbsp;<A
        href="http://www.R-project.org/posting-">http://www.R-project.org/posting-</A>
        &nbsp;</DIV>
        <DIV>&gt; &nbsp;guide.html</DIV>
        <DIV>&nbsp;</DIV></BLOCKQUOTE></TD></TR></TBODY></TABLE></DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;</DIV></BODY></HTML>

--=====003_Dragon527446281311_=====--

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: Weibull distribution

Valentin Dimitrov-2
In reply to this post by Leaf Sun
Dear Leaf,

I modified your code as follows:

gamma.fun <- function(mu,sd,start=100)    
 {
f.fn <- function(alpha)
{abs(sd^2-mu^2/(gamma(1+1/alpha))^2*(gamma(1+2/alpha)-(gamma(1+1/alpha))^2))}
alpha <- optim(start, f.fn)
beta <- mu/gamma(1+1/alpha$par)
return(list=c(a=alpha$par,b=beta));
 }

Now it works properly.

First, I added an abs(). You tried to solve an
equation by means of the R-function optim(), which
finds a minimum. That's why you can find the solution
of f(x)=a through minimization of abs(f(x)-a).
Second, I deleted the optim-method BFGS from the
optim() function, because it is not appropriate in
this case. BFGS seeks to make the gradient (or here
the first derivative) zero and in your case f(x)
converges to a constant for big x, which means f'(x)
is approximately 0 for big x, which is why BFGS stops
almost immediately after the start value. The default
method of optim() ( Nelder and Mead ) is more
appropriate, since it does not need the first
derivative and works only with function values.

Best regards,
Valentin


--- Leaf Sun <[hidden email]> wrote:

> Hi all,
>
> By its definition, the mean and variance of two-par.
> Weibull distribution are:
>
>  
>
>  
>
>  (www.wikipedia.org)
>
>
> I was wondering, if given mean and sd. could we
> parameterize the distribution? I tried this in R.
>
> gamma.fun <- function(mu,sd,start=100)    
> {
> f.fn <- function(alpha)
>
sd^2-mu^2/(gamma(1+1/alpha))^2*(gamma(1+2/alpha)-(gamma(1+1/alpha))^2)

> alpha <- optim(start, f.fn,method='BFGS')
> beta <- mu/gamma(1+1/alpha$par)
> return(list=c(a=alpha$par,b=beta));
> }
>
>
> But the problems come up here:
>
> 1)  the return values of a and b are only related to
> the input mean, and nothing to do with the sd. For
> instance, when I apply a mean mu = 3 whatever I use
> sd=2, sd=4, the function returned the same scale and
> shape values.
>
> > gamma.fun(3,4,10);
>        a        b
> 5.112554 3.263178
>
> > gamma.fun(3,2,10);
>        a        b
> 5.112554 3.263178
>
> 2) the start value determines the results: if I
> apply mean = 3, and sd=2, with a start of 10, it
> would return alpha close to 10, if I use a start =
> 100, it would return alpha close to 100.
>
> > gamma.fun(3,2,10);
>        a        b
> 5.112554 3.263178
>
> > gamma.fun(3,2,100);
>         a         b
> 99.999971  3.017120
>
> Since I am not a statistician, I guess there must be
> some theoretical reasons wrong with this question.
> So I am looking forward to some correction and
> advice to solve these. Thanks a lot in advance!
>
> Leaf
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
>

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: Weibull distribution

Thomas Lumley
On Fri, 21 Jul 2006, Valentin Dimitrov wrote:

> Dear Leaf,
>
> I modified your code as follows:
>
> gamma.fun <- function(mu,sd,start=100)
> {
> f.fn <- function(alpha)
> {abs(sd^2-mu^2/(gamma(1+1/alpha))^2*(gamma(1+2/alpha)-(gamma(1+1/alpha))^2))}
> alpha <- optim(start, f.fn)
> beta <- mu/gamma(1+1/alpha$par)
> return(list=c(a=alpha$par,b=beta));
> }
>
> Now it works properly.
>
> First, I added an abs(). You tried to solve an
> equation by means of the R-function optim(), which
> finds a minimum. That's why you can find the solution
> of f(x)=a through minimization of abs(f(x)-a).
> Second, I deleted the optim-method BFGS from the
> optim() function, because it is not appropriate in
> this case.

optim() is not appropriate at all in this case -- its help page says to
use optimize() for one-dimensional problems.

In fact, in one dimension there isn't any need to resort to optimization
when you really want root-finding, and uniroot() is more appropriate than
optimize().


  -thomas

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: Weibull distribution

Leaf Sun
Thanks for the suggestion! I switched to optimize(), al <- optimize(f.fn, lower = 0.1, upper =100,tol=0.001);
the warnings were gone and it works stably.
But when I tried  al <- uniroot(f.fn, lower = 0.1, upper =100,tol=0.001);
error occured: f() values at end points not of opposite sign. The error seems to me like there is no root found within the interval. I was not able to solve this problem.

Thanks!

Leaf





----- Original Message -----

From: Thomas Lumley,  [hidden email]
Sent: 2006-07-21,  09:35:11
To: Valentin Dimitrov, [hidden email]
Subject:  Re: [R] Weibull distribution
 
On  Fri,  21  Jul  2006,  Valentin  Dimitrov  wrote:

>  Dear  Leaf,
>
>  I  modified  your  code  as  follows:
>
>  gamma.fun   <-  function(mu,sd,start=100)
>  {
>  f.fn   <-  function(alpha)
>  {abs(sd^2-mu^2/(gamma(1+1/alpha))^2*(gamma(1+2/alpha)-(gamma(1+1/alpha))^2))}
>  alpha   <-  optim(start,  f.fn)
>  beta   <-  mu/gamma(1+1/alpha$par)
>  return(list=c(a=alpha$par,b=beta));
>  }
>
>  Now  it  works  properly.
>
>  First,  I  added  an  abs().  You  tried  to  solve  an
>  equation  by  means  of  the  R-function  optim(),  which
>  finds  a  minimum.  That's  why  you  can  find  the  solution
>  of  f(x)=a  through  minimization  of  abs(f(x)-a).
>  Second,  I  deleted  the  optim-method  BFGS  from  the
>  optim()  function,  because  it  is  not  appropriate  in
>  this  case.

optim()  is  not  appropriate  at  all  in  this  case  --  its  help  page  says  to  
use  optimize()  for  one-dimensional  problems.

In  fact,  in  one  dimension  there  isn't  any  need  to  resort  to  optimization  
when  you  really  want  root-finding,  and  uniroot()  is  more  appropriate  than  
optimize().


  -thomas

        [[alternative HTML version deleted]]

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: Weibull distribution

Valentin Dimitrov-2
It seems to me that not for all values of mu and sd
there is a Weibull distribution with mean=mu and
variance=sd^2.
the programm with optimize(f.fn) finds always a
solution, but this is not necessarily what we need,
because the minimum of (abs(f(x)) is not always 0.
Suppose f(x)=2+x^2, then optimize(x) finds x=0, but
x=0 is not a root of f(x)=0.
That's why I agree with Thomas Lumley, that uniroot()
could be more appropriate than optim and optimize.

Best regards,
Valentin

--- Leaf Sun <[hidden email]> wrote:

> Thanks for the suggestion! I switched to optimize(),
> al <- optimize(f.fn, lower = 0.1, upper
> =100,tol=0.001);
> the warnings were gone and it works stably.
> But when I tried  al <- uniroot(f.fn, lower = 0.1,
> upper =100,tol=0.001);
> error occured: f() values at end points not of
> opposite sign. The error seems to me like there is
> no root found within the interval. I was not able to
> solve this problem.
>
> Thanks!
>
> Leaf
>
>
>
>
>
> ----- Original Message -----
>
> From: Thomas Lumley,  [hidden email]
> Sent: 2006-07-21,  09:35:11
> To: Valentin Dimitrov, [hidden email]
> Subject:  Re: [R] Weibull distribution
>  
> On  Fri,  21  Jul  2006,  Valentin  Dimitrov  wrote:
>
> >  Dear  Leaf,
> >
> >  I  modified  your  code  as  follows:
> >
> >  gamma.fun   <-  function(mu,sd,start=100)
> >  {
> >  f.fn   <-  function(alpha)
> >
>
{abs(sd^2-mu^2/(gamma(1+1/alpha))^2*(gamma(1+2/alpha)-(gamma(1+1/alpha))^2))}

> >  alpha   <-  optim(start,  f.fn)
> >  beta   <-  mu/gamma(1+1/alpha$par)
> >  return(list=c(a=alpha$par,b=beta));
> >  }
> >
> >  Now  it  works  properly.
> >
> >  First,  I  added  an  abs().  You  tried  to
> solve  an
> >  equation  by  means  of  the  R-function
> optim(),  which
> >  finds  a  minimum.  That's  why  you  can  find
> the  solution
> >  of  f(x)=a  through  minimization  of
> abs(f(x)-a).
> >  Second,  I  deleted  the  optim-method  BFGS
> from  the
> >  optim()  function,  because  it  is  not
> appropriate  in
> >  this  case.
>
> optim()  is  not  appropriate  at  all  in  this
> case  --  its  help  page  says  to  
> use  optimize()  for  one-dimensional  problems.
>
> In  fact,  in  one  dimension  there  isn't  any
> need  to  resort  to  optimization  
> when  you  really  want  root-finding,  and
> uniroot()  is  more  appropriate  than  
> optimize().
>
>
>   -thomas
>

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: Weibull distribution

Göran Broström
In reply to this post by Thomas Lumley
On 7/21/06, Thomas Lumley <[hidden email]> wrote:

> On Fri, 21 Jul 2006, Valentin Dimitrov wrote:
>
> > Dear Leaf,
> >
> > I modified your code as follows:
> >
> > gamma.fun <- function(mu,sd,start=100)
> > {
> > f.fn <- function(alpha)
> > {abs(sd^2-mu^2/(gamma(1+1/alpha))^2*(gamma(1+2/alpha)-(gamma(1+1/alpha))^2))}
> > alpha <- optim(start, f.fn)
> > beta <- mu/gamma(1+1/alpha$par)
> > return(list=c(a=alpha$par,b=beta));
> > }
> >
> > Now it works properly.
> >
> > First, I added an abs(). You tried to solve an
> > equation by means of the R-function optim(), which
> > finds a minimum. That's why you can find the solution
> > of f(x)=a through minimization of abs(f(x)-a).
> > Second, I deleted the optim-method BFGS from the
> > optim() function, because it is not appropriate in
> > this case.
>
> optim() is not appropriate at all in this case -- its help page says to
> use optimize() for one-dimensional problems.

Just to clarify: The help page says:

'optim' will work with one-dimensional 'par's, but the default
method does not work well (and will warn).  Use 'optimize' instead.

In other words, if you for instance use the 'BFGS' method, optim is
perfectly OK for one-dimensional problems.
>
> In fact, in one dimension there isn't any need to resort to optimization
> when you really want root-finding, and uniroot() is more appropriate than
> optimize().

One reason to use optim instead of uniroot or optimize is that you
need not specify a finite interval that covers the solution.

Göran

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.