Il modello è troppo parametrizzato: non è necessario $ \ beta_1 $ , che può essere impostato su qualcosa di conveniente, come 1.

Una cosa a cui ho pensato è stata di adattarmi in modo iterativo. Inizia con qualche ipotesi su $ w $ e $ \ beta_2 $ . Quindi calcola $ Z = (\ sum_i \ hat {w} _iX_i) ^ 2 $ e adatta il modello lineare

Y ~ X1 + X2 + ... + X_k + Z

I coefficienti dei $ X $ sono i nuovi $ \ hat {w} _i $ e il coefficiente di $ Z $ è $ \ hat \ beta_2 $ . E poi ricalcola Z , itera e spera che converga. Purtroppo non è così.

Ma se $ k $ non è troppo grande, è facile calcolare semplicemente la somma residua dei quadrati in funzione dei parametri ed eseguirla attraverso un generale ottimizzatore di scopo. In R userei minqa :: newuoa , ma ci sono molte alternative.

  > X<-matrix (rnorm (50 * 100), ncol = 5)
> w<-1: 5
> Y<- (X% *% w) + 2 * (X% *% w) ^ 2 + rnorm (100)
>
>
> rss<-function (theta) {
+ beta2<-theta [1]
+ w<-theta [-1]
+ mu<- (X% *% w) + beta2 * (X% *% w) ^ 2
+ somma ((Y-mu) ^ 2)
+}
>
> minqa :: newuoa (par = rep (1,6), rss)
stime dei parametri: 1.99478699135839, 1.00032043499982, 2.00140284432351, 3.00312315850919, 4.00284240744153, 5.00537517104468
obiettivo: 1047.51402563294
numero di valutazioni di funzione: 1689
 

Quindi utilizza il bootstrap per ottenere stime di errore standard.

Con $ k = 50 $ non funziona (senza ottimizzazione - sono sicuro che funzionerebbe se le impostazioni predefinite dell'ottimizzatore fossero cambiate o i valori iniziali erano migliori)

Se scrivi l'espressione, ottieni un polinomio in termini di $ X_1, X_2, .., X_k $ , comprese le loro interazioni, dove il nuovo " coefficienti "sono tutti funzione di $ \ beta $ se $ w $ se due. Per k = 2, ottieni un polinomio che ha 5 coefficienti (o 6 compresa l'intercetta) con 4 incognite:

$$ \ begin {align *} Y & = \ beta_0 + (\ beta_1w_1) X_1 + (\ beta_1w_2) X_2 + (\ beta_2w_1 ^ 2) X_1 ^ 2 + (\ beta_2 w_2 ^ 2) X_2 ^ 2 + (2 \ beta_2 w_1w_2) X_1X_2 + \ varepsilon \\ & = \ alpha_0 + \ alpha_1X_1 + \ alpha_2X_2 + \ alpha_3X_1 ^ 2 + \ alpha_4X_2 ^ 2 + \ alpha_5X_1X_2 + \ varepsilon

Se soddisfi questa regressione, otterrai i nuovi coefficienti $ \ alpha $ , che ti danno un sistema di equazioni non lineari:

$$ \ begin {align *} \ alpha_0 & = \ beta_0 \\ \ alpha_1 & = \ beta_1w_1 \\ \ alpha_2 & = \ beta_1w_2 \\ \ alpha_3 & = \ beta_2w_1 ^ 2 \\ \ alpha_4 & = \ beta_2 w_2 ^ 2 \\ \ alpha_5 & = 2 \ beta_2 w_1w_2 \ end {align *} $$

In linea di principio, quel sistema di equazioni dovrebbe essere risolvibile numericamente, almeno a volte. Dovrebbe rimanere risolvibile con $ k>3 $ poiché non hai la maledizione della dimensionalità poiché ogni nuova variabile aggiunge solo un parametro ma più nuove equazioni che aiutano a fissarlo.

Ecco un esempio di simulazione giocattolo $ k = 2 $ utilizzando Stata in cui ignoro l'equazione dell'intercetta poiché è banale:

 . chiaro

. imposta obs 1000
il numero di osservazioni (_N) era 0, ora 1.000

. impostare seme 10011979

. gen b0 = 1

. gen b1 = 2

. gen b2 = 3

. gen w1 = 4

. gen w2 = 5

. gen x1 = normale (0,1)

. gen x2 = rnormale (10,2)

. gen eps = rnormal ()

. gen y = b0 + b1 * (w1 * x1 + w2 * x2) + b2 * (w1 * x1 + w2 * x2) ^ 2 + eps
. reg y (c.x1 c.x2) ## (c.x1 c.x2)

      Fonte | SS df MS Numero di oss = 1.000
------------- + ---------------------------------- F ( 5, 994) > 99999.00
       Modello | 1.1237e + 10 5 2.2475e + 09 Prob > F = 0,0000
    Residuo | 1052.11816 994 1.05846897 R quadrato = 1.0000
------------- + ---------------------------------- Adj R -squared = 1.0000
       Totale | 1.1237e + 10999 11248523.6 Root MSE = 1.0288

-------------------------------------------------- ----------------------------
           y | Coef. Std. Err. t P> | t | [95% Conf. Intervallo]
------------- + ------------------------------------ ----------------------------
          x1 | 8.082131 .1573906 51.35 0.000 7.773275 8.390987
          x2 | 9.852645 .110114 89.48 0.000 9.636562 10.06873
             |
   c.x1 # c.x1 | 47.9813 .0233895 2051.40 0.000 47.9354 48.0272
             |
   c.x1 # c.x2 | 119.9907 .0153233 7830.59 0.000 119.9606 120.0208
             |
   c.x2 # c.x2 | 75.00664 .0053927 1.4e + 04 0.000 74.99605 75.01722
             |
       _cons | 1.77947 .5532575 3.22 0.001 .693783 2.865156
-------------------------------------------------- ----------------------------

.
. mata chiaro

. mata:
------------------------------------------------- mata (digita end per uscire) -------------------------------------------- -------------------------------------------------- -------------------------------------------------
: void mysolver (todo, p, lnf, S, H)
> {
> b1 = p [1]
> b2 = p [2]
> w1 = p [3]
> w2 = p [4]
> lnf = (b1 * w1 - 8.082131) ^ 2 \
> (b1 * w2 - 9.852645) ^ 2 \
> (b2 * w1 ^ 2 - 47,9813) ^ 2 \
> (b2 * w2 ^ 2 - 75.00664) ^ 2 \
> (2 * b2 * w1 * w2 - 119.9907) ^ 2
>}
nota: argomento da fare inutilizzato
nota: argomento S inutilizzato
nota: argomento H inutilizzato

:
: S = ottimizzare_init ()

: optim_init_evaluator (S, &mysolver ())

: optim_init_evaluatortype (S, "v0")

: optim_init_params (S, (1,1,1,1))

: optim_init_which (S, "min")

: optim_init_tracelevel (S, "nessuno")

: optim_init_conv_ptol (S, 1e-16)

: optim_init_conv_vtol (S, 1e-16)

: p = ottimizza (S)

: p
                 1 2 3 4
    + ------------------------------------------------- -------- +
  1 | 2.1561597 3.521534782 3.691630188 4.614939185 |
    + ------------------------------------------------- -------- +

: fine
-------------------------------------------------- -------------------------------------------------- -------------------------------------------------- --------------------------------------------
 

La soluzione non è molto buona (a meno che non strizzi gli occhi e arrotondi all'intero più vicino), poiché $ p = (2,3,4,5) $ in la simulazione. Probabilmente sto facendo qualcosa di sbagliato quando risolvo numericamente le equazioni. Ma anche l'intercetta è piuttosto fuori luogo con $ b_0 = 1.77947 \ ne 1 $ .

Code:

  cls
chiaro
imposta obs 1000
impostare seme 10011979
gen b0 = 1
gen b1 = 2
gen b2 = 3
gen w1 = 4
gen w2 = 5
gen x1 = normale (0,1)
gen x2 = rnormale (10,2)
gen eps = rnormal ()
gen y = b0 + b1 * (w1 * x1 + w2 * x2) + b2 * (w1 * x1 + w2 * x2) ^ 2 + eps
reg y (c.x1 c.x2) ## (c.x1 c.x2)

mata chiaro
mata:
void mysolver (todo, p, lnf, S, H)
         {
                 b1 = p [1]
                 b2 = p [2]
                 w1 = p [3]
                 w2 = p [4]
                 lnf = (b1 * w1 - 8.082131) ^ 2 \
(b1 * w2 - 9,852645) ^ 2 \
                       (b2 * w1 ^ 2 - 47,9813) ^ 2 \
                       (b2 * w2 ^ 2 - 75.00664) ^ 2 \
                       (2 * b2 * w1 * w2 - 119,9907) ^ 2
        }

S = ottimizzare_init ()
ottimizzare_init_evaluator (S, &mysolver ())
optim_init_evaluatortype (S, "v0")
ottimizzare_init_params (S, (1,1,1,1))
optim_init_which (S, "min")
ottimizzare_init_tracelevel (S, "nessuno")
ottimizzare_init_conv_ptol (S, 1e-16)
ottimizzare_init_conv_vtol (S, 1e-16)
p = ottimizza (S)
p
fine