Maximum Likelihood – How to Use R to Plot $f_n(\theta)$ from Samples Effectively [Closed]

Question

This question appears to be off-topic because it focuses on programming, debugging, or performing routine operations, or it asks about obtaining datasets. You could try the support links we maintain or the Open Data site instead.

---

If the question is actually a statistical topic disguised as a coding question, then OP should edit the question to clarify this. After the statistical content has been clarified, the question is eligible for reopening.

Closed 2 years ago.

Improve this question

For iid random samples $X_1,\dots, X_n$ with mean value $E X_1=\theta$, take $X_{(1)}=\min \{X_1,\dots, X_n\}$ and $X_{(n)}=\max \{X_1,\dots, X_n\}$. If $\lambda$ solves equation
$$
\sum_{i=1}^n \frac{X_i-\theta}{1+\lambda(X_i-\theta)}=0
$$
for $\lambda\in [\frac{1/n-1}{X_{(n)-\theta}}, \frac{1/n-1}{X_{(1)-\theta}}]$.

Let
$$
f_n(\theta)=-\sum_{i=1}^n\log(1+\lambda(X_i-\theta))
$$

How to use R to plot $f_n(\theta)$ with $\theta\in [1,12]$ based on following code?

# I first get 100 samples 
set.seed(201111)
x=rlnorm(100,0,2)

#theta is the mean of x
theta=mean(x)
n=100
xmax=max(x)
xmin=min(x)
# I define function h, and use unifoot function to find lambda
h=function(lam, n)
{
 sum((x-theta)/(1+lam*(x-theta)))
}
uniroot(h, c((1/n-1)/(xmax-theta),(1/n-1)/(xmin-theta)),n=n)$root

But how to plot `f_n(theta)' like ?

Based on the code in answer, I got the following point plot.

Answer 1

You can create a sequence that goes from 1 to 12 at whatever granularity you desire, use a for loop to compute 'f_n(theta)' at each point in the sequence, storing the results in a vector. Then at the end plot that vector, and it will give you a plot of 'f_n(theta)' over [1,12]. You can create a sequence with the seq() function.

theta_seq = seq(from = 1, to = 12, by = 0.01)

I'd also ask though, is $\theta$ supposed to be $E(X_1)$ or $\bar{X}_n$? If the former, I'd interpret your example as needing to have

theta = 0

rather than

theta = mean(x)

You could ultimately get something like this

set.seed(201111)
theta_seq = seq(from = 1, to = 12, by = 0.01)
f_n = rep(NA, length(theta_seq))

# I define function h, and use unifoot function to find lambda
h=function(lam, n)
{
 sum((x-theta)/(1+lam*(x-theta)))
}

f_n = rep(NA, length(theta_seq))
i = 1
x = rnorm(100, 0, 2)
n=100
xmax=max(x)
xmin=min(x)

for (theta in theta_seq)
{
lambda = uniroot(h, c((1/n-1)/(xmax-theta),(1/n-1)/(xmin-theta)),n=n)$root

f_n[i] = -sum(log(1+lambda*(x-theta)))
i = i+1
}
plot(theta_seq, f_n, type = "l")

I get a plot that looks like this