1
$\begingroup$

Equations :

I'm trying to find the best way of estimating the parameters of the following equations :

  • Knowing the value of $n$ (in my case $n=27$),

  • Knowing the values of constants $(a_i)_{i \in \{1,...,n\}}$, $(b_j)_{j \in \{1,...,n\}}$ and $(c_k)_{i \in \{1,...,n\}}$, both in $\mathbb{R}$

I want to solve the following system of $81$ equations for $(x_i)_{i \in \{1,...,n\}}$, $(y_j)_{j \in \{1,...,n\}}$ and $(z_k)_{i \in \{1,...,n\}}$, both in $\mathbb{R}$

$$\forall i \in \{1,...,n\}, a_i = \sum_{j = 1}^{n+1-i} \big{[} \sum_{k = 1}^{n+2-i-j} x_i y_j z_k \big{]}$$

$$\forall j \in \{1,...,n\}, b_j = \sum_{i = 1}^{n+1-j} \big{[} \sum_{k = 1}^{n+2-j-i} x_i y_j z_k \big{]}$$

$$\forall k \in \{1,...,n\}, c_j = \sum_{i = 1}^{n+1-k} \big{[} \sum_{j = 1}^{n+2-k-i} x_i y_j z_k \big{]}$$

For an algorythmic way of solving thoose, see remark below

Thoose equations correspond to the score of my model. I'm of course able to derive them.

Context :

In a quasipoisson GLM with 3 sets of qualitative variables (respectively indexed by $i$,$j$ and $k$ ), thoose equations correpond to the score of the model if you choose $x_i = e^{a_i}$ and $y$ and $z$ accordingly to reparametrise. This reparametrisation allow for negative parameters. The following equations, for $a$, $b$ and $c$ each in $\mathbb{R}_{+}^{n}$, would then be a proper score for a quasi-poisson glm :

$$\forall i \in \{1,...,n\}, a_i = \sum_{j = 1}^{n+1-i} \big{[} \sum_{k = 1}^{n+2-i-j} e^{\alpha_i +\beta_j +\gamma_k} \big{]}$$

$$\forall j \in \{1,...,n\}, b_j = \sum_{i = 1}^{n+1-j} \big{[} \sum_{k = 1}^{n+2-j-i} e^{\alpha_i +\beta_j +\gamma_k} \big{]}$$

$$\forall k \in \{1,...,n\}, c_j = \sum_{i = 1}^{n+1-k} \big{[} \sum_{j = 1}^{n+2-k-i} e^{\alpha_i +\beta_j +\gamma_k} \big{]}$$

In the case of my data, i have and i want to fit the negatives values in this context. Of course thoose equations are not anymore MLE scores equations of a quasi-glm log-poisson model !

Remark :

One could notice that the solutions of such equations are not unique up to 2 constants : if you set $x_2 = fg x$, $y_2 = f y$ and $z_2 = g z, \forall(f,g) \in \mathbb{R}^2$, the equations are still the same.

So 2 more constraint could be added to the problem. The 2 logical constraint taking account of (complicated) context could be :

$$\sum_{j=1}^{n} y_j = 1$$ $$\sum_{k=1}^{n} z_j = 1$$

But since a re-parametrisation is easy to do, you could choose anything else if you prefer.

What was tried :

I'm working with R.

Derivation of the equations

The jacobian was easy to calculate since the derivation is quite straightforward.

R-base::optim

I tried to do a least-sum-of-square evaluation, and minimising the sum-of-square with optim.

  1. It worked well, and gave me values that i was expecing for.

  2. It took more than 8 hours to get to the optimum.

Since i'm gonna have to solve thoose kind of system $N = 10000$ times ( in a bootstrap context), this type of calculation time is'nt acceptable. However, i do not know how to speed up things.

nleqslv::nleqslv

Giving the jacobian and my equations to $nleqslv$ allowed me to check ( via the chkJac option) that my implementation was right. However, the time of computation was still very high... The algorythma did converge to the right optimum too.

R-code implementing the equations.

If you want to try, i give here the R code implementing the equations, and a sample of possible constants.

# le score du problème : 
.score <- function(par,obs=observed_score){
  score = vector(length=81)

  for(i in 1:27){
    for(j in 28:(54-i+1)){
      for(k in 55:(81-i-(j-27)+2)){
        score[i] <- score[i] + par[i] * par[j] * par[k]
        score[j] <- score[j] + par[i] * par[j] * par[k]
        score[k] <- score[k] + par[i] * par[j] * par[k]
      }
    }
  }
  return(score - obs)
}


.jacobian_col <- function(par,der){
  score = rep(0,81)
  for(i in 1:27){
    for(j in 28:(54-i+1)){
      for(k in 55:(81-i-(j-27)+2)){
        score[i] <- score[i] + prod(par[c(i,j,k)[c(i,j,k)!=der]])*(der%in%c(i,j,k))
        score[j] <- score[j] + prod(par[c(i,j,k)[c(i,j,k)!=der]])*(der%in%c(i,j,k))
        score[k] <- score[k] + prod(par[c(i,j,k)[c(i,j,k)!=der]])*(der%in%c(i,j,k))
      }
    }
  }
  return(score)
}

.jacobian <- function(par){
  do.call(cbind,map(seq_along(par),~.jacobian_col(par,.x)))
}

structure(c(11570.1643291449, 10421.1383991299, 15354.8942350771, 
14672.5980606024, 25625.9411999004, 31961.0161884923, 43378.9129173754, 
67814.381008179, 60714.7908481998, 39529.4222899274, 27887.1320742589, 
42230.2912306984, 37897.3415405366, 42871.5838043031, 62596.8779588982, 
61718.6045656741, 53139.9646059978, 38056.8547769802, 33021.286638975, 
37135.4947146619, 33909.3051389397, 23701.8900488388, 14757.3839741324, 
11899.1805782697, 6873.14393317911, 2126.29153481316, 679.052493645299, 
25756.7273691109, 99813.2137201322, 109786.936358354, 118185.697766755, 
104666.918792329, 89897.5139872207, 72800.2957355428, 60348.8765301957, 
47393.6293273173, 47846.1322946524, 42593.2294177613, 20983.0205448836, 
7233.54539327032, 1686.16523097222, 1527.64863388384, 335.295487558557, 
193.158288739573, 0.414306553170475, 456.51193749324, 4.82376149907932, 
0.0173205080756888, 0, 0, 0, 25.9824941643407, 9.1843899327228, 
0, 428326.791718838, 132378.733653794, 73766.9437996086, 60834.4178297494, 
53800.0175068735, 27983.4955869634, 28557.5423618357, 19401.7642630057, 
14580.3733525847, 3311.19528649966, 3652.29368119933, -4709.2016841288, 
4982.91729219773, -518.174279921813, -305.57174632936, 6000.83975861469, 
-66.6821374380533, -629.658422971308, -0.859253085127191, -31.3862501968426, 
58.4934169120623, 34.4323386700816, 137.739902869745, 0.195842984811725, 
0.0173205080756888, -1.73205080756888, 0), .Dim = c(27L, 3L), .Dimnames = list(
NULL, c("margin.sum.i", "margin.sum.j", "margin.sum.k")))

My question / TL;DR

I wanna solve numericaly the equation from the first point of this post, in the fastest time possible. I provided in the code the constants you'll need and the jacobian of the eqations;

I'm open to any suggestions.

$\endgroup$
2
  • $\begingroup$ Not sure how to solve your system of equations, but if you want to estimate a Maximum Likelihood problem without an analytical likelihood function, you could look at adversarial techniques like a Generative Adversarial Network (GAN). $\endgroup$ – www3 Mar 1 '18 at 20:02
  • $\begingroup$ It is not anymore a likelyhood problem since the parameter ranges are outside what the likelyhood or quasi-likelyhood would ask for. So thoose kind of methods wont work.. $\endgroup$ – Pilou l'Nrv Mar 2 '18 at 8:04

Your Answer

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

Browse other questions tagged or ask your own question.