find MLE od theta from cauchi distrubtion with Newton Rapson method

General
1

I have this sample data -

> dput(my_vec)
c(7.16523478153752, 5.66659652818595, 4.47575534893755, 4.84970857977856, 
15.2276296414708, -0.573093658844655, 4.97980673868322, 2.73969325233614, 
5.14683035133365, 10.1221488713611, 9.01656611721311, 65.711819422978, 
5.52057043354834, 6.30674880627702, 8.67771771267678, 5.2528503049587, 
3.50395623925858, 4.24774012371174, 11.4137624410312, -48.1722033880239, 
-0.376400642113356, 5.76475359419462, -27.353313803102, 4.09682042042852, 
5.03373747558625, 3.8261660769812, 4.43580895525249, 4.22242932760446, 
4.44905425775097, 4.98475525084258, 3.6416524979406, 3.81767927422987, 
-93.2141354888589, 5.01103555428068, 5.38206564752185, 3.0296536606134, 
86.676038167071, 3.76536864898752, 7.26590572567999, 4.63759338498908, 
-17.2259677581435, 4.90903933854569, 4.43042810097151, 5.41519101700693, 
7.01553803666926, 5.05370712380939, 5.22896180532498, 3.92429687923716, 
5.46452912130986, 4.50524864464907, 6.13119216629816, 6.931011365041, 
5.70039361940988, 6.12470771804837, 6.66119827017415, -4.26865103703534, 
4.77160581324527, 7.91297525486072, 7.04882594451997, -98.2224152538262, 
6.22663920777969, 5.77535246011091, -9.91868097075743, 7.77803716672203, 
-10.1297033914588, 4.56699263921898, 8.56120643036614, 2.28239965661689, 
5.6212927060249, 5.42419784358529, 5.3654428914847, 3.89730116338776, 
3.93504254066386, 4.38168766024675, 3.00038017933155, 4.81884527713388, 
4.45257297902316, -3.50516232920359, 6.07365636649988, 4.26195094225287, 
4.73753258557777, 0.807607628402986, 3.93740907315333, 3.08283749617447, 
3.77379774203008, 2.56562208889432, -19.6532587812275, 8.00379422844706, 
5.27350155029373, 5.17570743606727, 6.49856446395129, -8.78462344722329, 
4.38775671122269, 4.39685118092286, 3.52571990926583, 7.13834590807858, 
0.724655622024244, 5.72807728660989, 6.58172179467414, 6.22426074825281
)

and, I need to estimate theta from our sample ( our data is the sample ) , with Newton Rapson method. theta parameter comes from a Cauchi distribution, when my initial value is the mean of the sample.

I wrote something like this -

> data<-read.table(file.choose(), header = TRUE, sep= "")
> my_vec<-as.numeric(as.character(unlist(data[[1]])))
> dput(my_vec)
> nr_method<-function(x,toler=.001){      #x is a vector here
>   start_value<-median(x)
>   n=length(x);
>   theta=start_value;
>   # Compute first deriviative of log likelihood
>   first_derivate<-2*sum((x-theta)/(1+(x-theta)^2))
>   # Continue Newton’s method until the first derivative
>   while(abs(first_derivate)>toler){
>     # Compute second derivative of log likelihood
>     second_derivate=2*sum(((x-theta)^2-1)/(1+(x-theta)^2)^2);
>     # Newton’s method update of estimate of theta
>     theta_new<-theta-first_derivate/second_derivate;
>     theta=theta_new;
>     # Compute first derivative of log likelihood
>     first_derivate=2*sum((x-theta)/(1+(x-theta)^2))
>   }
>   list(theta_hat=theta);
> }
> main_fun <- function (theta, x)
> {
>   sum (-log(1+(x-theta)^2))
> }
> th<-seq(-15,15,by=0.1)
> ll<-rep(0,length(th))
> for (i in 1:length(th)) ll[i] <- main_fun(th[i],x)
> plot(th,ll,type="l")
> try1<-nr_method(x,1e-5)

but I don't know ho to contunie from here.
Finanlly, I neeed to plot the likelihood, and make sure I found the max.

any help?

2

I see a few times in your code where you have x but need to have my_vec. I'm not sure what else you need to fix than that:

my_vec <- c(7.16523478153752, 5.66659652818595, 4.47575534893755, 4.84970857977856, 
            15.2276296414708, -0.573093658844655, 4.97980673868322, 2.73969325233614, 
            5.14683035133365, 10.1221488713611, 9.01656611721311, 65.711819422978, 
            5.52057043354834, 6.30674880627702, 8.67771771267678, 5.2528503049587, 
            3.50395623925858, 4.24774012371174, 11.4137624410312, -48.1722033880239, 
            -0.376400642113356, 5.76475359419462, -27.353313803102, 4.09682042042852, 
            5.03373747558625, 3.8261660769812, 4.43580895525249, 4.22242932760446, 
            4.44905425775097, 4.98475525084258, 3.6416524979406, 3.81767927422987, 
            -93.2141354888589, 5.01103555428068, 5.38206564752185, 3.0296536606134, 
            86.676038167071, 3.76536864898752, 7.26590572567999, 4.63759338498908, 
            -17.2259677581435, 4.90903933854569, 4.43042810097151, 5.41519101700693, 
            7.01553803666926, 5.05370712380939, 5.22896180532498, 3.92429687923716, 
            5.46452912130986, 4.50524864464907, 6.13119216629816, 6.931011365041, 
            5.70039361940988, 6.12470771804837, 6.66119827017415, -4.26865103703534, 
            4.77160581324527, 7.91297525486072, 7.04882594451997, -98.2224152538262, 
            6.22663920777969, 5.77535246011091, -9.91868097075743, 7.77803716672203, 
            -10.1297033914588, 4.56699263921898, 8.56120643036614, 2.28239965661689, 
            5.6212927060249, 5.42419784358529, 5.3654428914847, 3.89730116338776, 
            3.93504254066386, 4.38168766024675, 3.00038017933155, 4.81884527713388, 
            4.45257297902316, -3.50516232920359, 6.07365636649988, 4.26195094225287, 
            4.73753258557777, 0.807607628402986, 3.93740907315333, 3.08283749617447, 
            3.77379774203008, 2.56562208889432, -19.6532587812275, 8.00379422844706, 
            5.27350155029373, 5.17570743606727, 6.49856446395129, -8.78462344722329, 
            4.38775671122269, 4.39685118092286, 3.52571990926583, 7.13834590807858, 
            0.724655622024244, 5.72807728660989, 6.58172179467414, 6.22426074825281
)
nr_method<-function(x,toler=.001){      #x is a vector here
  start_value<-median(x)
  n=length(x);
  theta=start_value;
  # Compute first deriviative of log likelihood
  first_derivate<-2*sum((x-theta)/(1+(x-theta)^2))
  # Continue Newton's method until the first derivative
  while(abs(first_derivate)>toler){
    # Compute second derivative of log likelihood
    second_derivate=2*sum(((x-theta)^2-1)/(1+(x-theta)^2)^2);
    # Newton's method update of estimate of theta
    theta_new<-theta-first_derivate/second_derivate;
    theta=theta_new;
    # Compute first derivative of log likelihood
    first_derivate=2*sum((x-theta)/(1+(x-theta)^2))
  }
  list(theta_hat=theta);
}
main_fun <- function (theta, x)
{
  sum (-log(1+(x-theta)^2))
}
th<-seq(-15,15,by=0.1)
ll<-rep(0,length(th))
for (i in 1:length(th)) ll[i] <- main_fun(th[i],my_vec)
plot(th,ll,type="l")


try1<-nr_method(my_vec,1e-5)
try1                
#> $theta_hat
#> [1] 4.926955

Created on 2020-12-21 by the reprex package (v0.3.0)

3

yes, but the plot shouldn't look like this. it should be look like normal distrubtuion and with xlimit = c ( 4.0 , 5.4 ) and on ylimit = c ( 0, 2.5 )

but it doesn't work.

our initial value is the median of the sample . # this is the first step
then, we we need to update theta(k).
we need to use newton - raphson method and to run this step to convergence - until the difference between Theta(k-1) to Theta(k) is a negligible difference. ( for example 10^-4) .
we need to return the last Theta(k) as Theta_hat MLE

4

when I changed the limit it doesn't looks good. so, maybe I have miss something

5

I'm not sure why this will look like the normal distribution. It is the Cauchy distribution. Please refer to the homework policy - this is as far as I feel comfortable helping.

6

I have written a code, I just don't understand what is the problem in my code. I just asked for someone to look on my code. that's all.

7

puttting aside the homework policy controversy,.. how would you answer StatSteph's seemingly reasonable query?

I'm not sure why this will look like the normal distribution. It is the Cauchy distribution.

8

it's need to be a little bit similar to normal.. I didn't say it's need to be normal. my plot is almost there.... but I need some little help...

9

This topic was automatically closed 21 days after the last reply. New replies are no longer allowed.

If you have a query related to it or one of the replies, start a new topic and refer back with a link.