ElectroYou - la comunità dei professionisti del mondo elettrico

da **RenzoDF** » 21 dic 2009, 16:12

Bravissimo =D>

Grazie ... ho trovato qualcos'altro da aggiungere ai compiti per le vacanze :mrgreen:

BTW dato che ci sei ... non potresti provare anche con gli ultimi 9 punti postati :mrgreen:

da **xyz** » 21 dic 2009, 18:39

RenzoDF ha scritto:non potresti provare anche con gli ultimi 9 punti postati

Trovare la funzione che approssima meglio i punti dati è difficile visto i pochi punti a disposizione.

La funziona con l'errore minimo che trovato è questa:

$f(x,y)=2.6947*\sin(-1.2053*x)*\cos(0.9129*y)+4.1078$

sicuramente esiste qualcosa di migliore.

Lo script per R:

Codice: Seleziona tutto: x=c(1,1,1,2,2,2,3,3,3) y=c(1,2,3,1,2,3,1,2,3) out=c(3.2,4.4,6.5,2.5,4.7,5.8,5.1,3.6,2.9) nls(out ~ a1*sin(w1*x)*cos(w2*y)+k, start = list(a1=1,w1=1,w2=1,k=-1))

risultato:

Codice: Seleziona tutto: Nonlinear regression model model: out ~ a1 * sin(w1 * x) * cos(w2 * y) + k data: parent.frame() a1 w1 w2 k 2.6947 -1.2053 0.9129 4.1078 residual sum-of-squares: 0.902

da **IsidoroKZ** » 21 dic 2009, 18:49

xyz ha scritto:La funziona con l'errore minimo che trovato è questa:

$f(x,y)=2.6947*\sin(-1.2053*x)*\cos(0.9129*y)+4.1078$

La funzione di tipo a1 * sin(w1 * x) * cos(w2 * y) + k l'hai scelta tu o l'ha presa il programma da un suo data base di funzioni, cercandone una che fosse ragionevole per quei dati?

da **xyz** » 21 dic 2009, 19:52

IsidoroKZ ha scritto:La funzione di tipo a1 * sin(w1 * x) * cos(w2 * y) + k l'hai scelta tu o l'ha presa il programma da un suo data base di funzioni, cercandone una che fosse ragionevole per quei dati?

Ho plottato i punti in uno spazio 3D e ho visto la loro distribuzione altalenante, ho provato con semplici funzioni lineari ma ottenevo degli errori di stima molto alti, poi sono passato con funzioni esponenziali e logaritmiche ma non davano risultati accertabili, quando sono poi passato a inventarmi funzioni trigonometriche è saltata fuori quella espressione.

In generale per trovare le funzioni che approssimano meglio i dati uno si basa molto sulla propria esperienza e servono molti punti per capire a priori quale strada prendere.

da **RenzoDF** » 21 dic 2009, 19:53

Io una soluzione la ho ... e non è trascendente ma polinomiale :!:

... e con residuo nullo :wink:

vi lascio indovinare

da **RenzoDF** » 22 dic 2009, 0:39

Mi sono perso il link dal quale avevo appreso il metodo nei suoi passi elementari ... e copiato i 9 punti di partenza :mrgreen:

(ma appena lo trovo lo posto) ...
comunque partendo anche da qua http://www.codecogs.com/reference/maths ... lation.php
e usando Vandermonde http://mathworld.wolfram.com/VandermondeMatrix.html ...
con l'aiuto di uno script di Scilab e inserendo le coordinate x,y,z dei 9 punti
(1, 1, 3.2), (1, 2, 4.4), (1, 3, 6.5), (2, 1, 2.5), (2, 2, 4.7), (2, 3, 5.8 ), (3, 1, 5.1), (3, 2, 3.6), (3, 3, 2.9).

: V1.png (15.65 KiB) Osservato 6122 volte

otterremo i coefficienti del polinomio interpolatore

: v2.png (12.8 KiB) Osservato 6120 volte

plottando la superficie interpolante i 9 punti

: v3.gif (15.15 KiB) Osservato 6120 volte

e ricavando infine l'equazione della stessa
$z=f(x,y)=0.\text{975}x^{2}\text{y}^{2}\text{-5}.\text{275}x^{2}\text{y+ 5}.\text{95x}^{2}\text{-3}.\text{925}xy^{2}\text{+19}.\text{825}xy\text{ -21}.\text{55}x\text{ + 3}.\text{4}y^{2}\text{-14}.\text{7}y\text{ +18}.\text{5}$

da **xyz** » 22 dic 2009, 10:43

Non mi è venuto in mente di provare le combinazioni incrociate tra x e y per i termini del polinomio #-o

lo script per R:

Codice: Seleziona tutto: x=c(1,1,1,2,2,2,3,3,3) y=c(1,2,3,1,2,3,1,2,3) out=c(3.2,4.4,6.5,2.5,4.7,5.8,5.1,3.6,2.9) nls(out ~ a1*x+a2*x^2+b1*y+b2*y^2+c11*x*y+c21*x^2*y+c12*x*y^2+c22*x^2*y^2+k, start = list(a1=1,a2=1,b1=1,b2=1,c11=1,c21=1,c12=1,c22=1,k=1), algorithm = "port")

risultato:

Codice: Seleziona tutto: Nonlinear regression model model: out ~ a1 * x + a2 * x^2 + b1 * y + b2 * y^2 + c11 * x * y + c21 * x^2 * y + c12 * x * y^2 + c22 * x^2 * y^2 + k data: parent.frame() a1 a2 b1 b2 c11 c21 c12 c22 k -21.550 5.950 -14.700 3.400 19.825 -5.275 -3.925 0.975 18.500 residual sum-of-squares: 4.385e-26

l'errore è molto basso.

Il polinomio è lo stesso:

$0.975\,x^2\,y^2-3.925\,x\,y^2+3.4\,y^2-5.275\,x^2\,y+19.825\,x\,y- 14.7\,y+5.95\,x^2-21.55\,x+18.5$

Anch'io uso Scilab (me lo compilo pure) se mi serve.

da **RenzoDF** » 23 dic 2009, 17:21

Grazie a xyz, un nuovo FreeTool per Electroportal :!:

Posto una versione ampliata del codice di xyz per ottenere anche il plottaggio 3D della superficie trovata

Codice: Seleziona tutto: #Dati iniziali x=c(1,1,1,2,2,2,3,3,3) y=c(1,2,3,1,2,3,1,2,3) out=c(3.2,4.4,6.5,2.5,4.7,5.8,5.1,3.6,2.9) #Ricerca parametri funzione renzodf<- nls(out ~ c10*x+c20*x^2+c01*y+c02*y^2+c11*x*y+c21*x^2*y+c12*x*y^2+c22*x^2*y^2+k, start = list(c10=1,c20=1,c01=1,c02=1,c11=1,c21=1,c12=1,c22=1,k=1), algorithm = "port") #Definizione funzione interpolante coefficienti<- coef(renzodf) valore<-function(x,y){cbind( x,x^2,y,y^2,x*y,x^2*y,x*y^2,x^2*y^2,1) %*% coefficienti} #Plotting 3D x<-y<-seq(0.5, 3.5, len = 20) z<-outer(x,y,valore) persp(x, y, z, main="Electroportal Interpolation", col="lightgreen", theta=30, phi=20, r=50, d=0.5, expand=0.5, ltheta=40, lphi=180,shade=0.75, cex.axis=0.5,cex.lab=0.7, xlab = "x", ylab = "y", zlab="z", ticktype="detailed",nticks=5)

con risultato

: R.png (17.91 KiB) Osservato 6033 volte

da **xyz** » 23 dic 2009, 23:29

La gestione dei plot 3D in R non è interattiva ma esiste il modulo rgl per usare le OpenGL in tempo reale:

http://cran.r-project.org/web/packages/rgl/index.html

Un esempio d'utilizzo con lo script precedente:

Codice: Seleziona tutto: #Plotting 3D (OpenGL) require(rgl) x<-y<-seq(0.5, 3.5, len = 20) z<-outer(x,y,valore) open3d() bg3d("white") view3d(theta=0, phi=-70) persp3d(x, y, z, xlab="x", ylab="y", zlab="z", col="lightgreen", back="line", smooth=FALSE, aspect=c(1, 1, 0.5)) #automatic spin #play3d(spin3d(axis=c(0,0,1), rpm=5), duration=10) #screen snapshot PNG #snapshot3d("screen.png", fmt="png") #screen snapshot EPS #rgl.postscript("screen.eps", fmt = "eps") #save animated GIF #movie3d(spin3d(axis=c(0,0,1), rpm=10), fps=10, duration=5, dir=".", movie="anim3d")

Oltre a ruotare i dati in tempo reale col mouse, permette di ruotare in automatico la viewport o salvare lo lo schermo nei formati PNG o EPS. E' possibile creare semplici animazioni GIF: