Trazar la distribución gaussiana en 3D

10

En la teoría de la probabilidad, la distribución normal (o gaussiana) es una distribución de probabilidad continua muy común. Las distribuciones normales son importantes en estadística y a menudo se usan en las ciencias naturales y sociales para representar variables aleatorias de valor real cuyas distribuciones no se conocen.

El reto

Su desafío es trazar la densidad de probabilidad de la distribución gaussiana en un plano tridimensional . Esta función se define como:

Dónde:




A = 1, σ x = σ y = σ

Reglas

  • Su programa debe tomar una entrada σ , la desviación estándar.
  • Su programa debe imprimir un diagrama 3D de la distribución gaussiana en la más alta calidad según lo permita su idioma / sistema.
  • Su programa no puede utilizar una distribución gaussiana directa o una densidad de probabilidad incorporada.
  • Su programa no tiene que finalizar.
  • Su trama puede ser en blanco y negro o color.
  • Su parcela debe tener líneas de cuadrícula en la parte inferior. Las líneas de cuadrícula en los lados (como se muestra en los ejemplos) son innecesarias.
  • Su diagrama no necesita tener números de línea al lado de las líneas de la cuadrícula.

Puntuación

Como es habitual en , ¡el envío con la menor cantidad de bytes gana! Puede que nunca "acepte" una respuesta usando el botón, a menos que sea increíblemente pequeño e intuitivo.

Salida de ejemplo

Su salida podría verse así:

5 5

O podría verse así:

6 6

Salidas más válidas . Salidas inválidas .

MD XF
fuente
Estaba confundido de que acaba de mostrar la función para el eje X. ¿Necesitamos tomar entradas / salidas separadas para los sigma X e Y y mu?
Scott Milner
Entonces, ¿debemos suponer que μ es igual a 0? ¿Y qué escala necesita para x e y? Si los rangos x e y se eligen muy pequeños en relación con σ, entonces el gráfico se verá esencialmente como una función constante.
Greg Martin
(Para la distribución bidimensional, creo que es más claro si usas | x-μ | ^ 2 en la definición en lugar de (x-μ) ^ 2.)
Greg Martin
@GregMartin Editado.
MD XF
2
Todavía no está claro ... ¿qué son x_o y y_o y θ?
Greg Martin

Respuestas:

7

Gnuplot 4, 64 62 61 60 47 bytes

(Atado con Mathematica ! WooHoo!)

se t pn;se is 80;sp exp(-(x**2+y**2)/(2*$0**2))

Guarde el código anterior en un archivo llamado A.gpe invoque con lo siguiente:

gnuplot -e 'call "A.gp" $1'>GnuPlot3D.png

donde se $1va a reemplazar con el valor de σ. Esto guardará un .pngarchivo llamado que GnuPlot3D.pngcontiene la salida deseada en el directorio de trabajo actual.

Tenga en cuenta que esto solo funciona con distribuciones de Gnuplot 4, ya que en Gnuplot 5 las $nreferencias a los argumentos fueron desaprobadas y reemplazadas por las más desafortunadamente más detalladas ARGn.

Salida de muestra con σ = 3:

Salida de muestra

Esta salida está bien según OP .


Gnuplot 4, solución alternativa, 60 bytes

Aquí hay una solución alternativa que es mucho más larga que la anterior, pero el resultado se ve mucho mejor en mi opinión.

se t pn;se is 80;se xyp 0;sp exp(-(x**2+y**2)/(2*$0**2))w pm

Esto todavía requiere Gnuplot 4 por la misma razón que la solución anterior.

Salida de muestra con σ = 3:

Salida de muestra n. ° 2

R. Kap
fuente
I am not sure if it molds to the specifications required¿Qué especificaciones crees que no cumple?
MD XF
@MDXF En primer lugar, no estoy seguro de si la transparencia del gráfico está bien. Sinceramente, no me gusta mucho, por eso no estaba seguro de si estaría bien aquí. En segundo lugar, el gráfico comienza a una unidad de altura desde la parte inferior de forma predeterminada, y tampoco estoy seguro de si eso está bien. En tercer lugar, debido a que el gráfico comienza una unidad de altura, no estoy seguro de que la desproporcionalidad del gráfico en comparación con los gráficos proporcionados en la publicación original sea correcta. Sin embargo, si todo está bien con usted, felizmente lo convertiré en la respuesta principal.
R. Kap
@MDXF De hecho, iba a publicarlo como la respuesta original, pero por estas razones elegí no publicarlo y respondí por respuesta actual.
R. Kap
@MDXF En realidad, puedo hacerlo aún más corto si esto está bien. Entiendo si no lo será, pero no está de más preguntar. Es la forma predeterminada de Gnuplotgraficar la densidad de probabilidad de la distribución gaussiana con un Sigma 2sin ninguna modificación del entorno.
R. Kap
@MDXF Creo que podría haber preguntado antes de publicar mi respuesta original, pero en ese momento estaba muy ansioso por publicar una respuesta.
R. Kap
14

C ++, 3477 3344 bytes

El recuento de bytes no incluye las nuevas líneas innecesarias.
MD XF jugó 133 bytes.

No hay forma de que C ++ pueda competir por esto, pero pensé que sería divertido escribir un procesador de software para el desafío. Arranqué y jugué algunos fragmentos de GLM para las matemáticas 3D y utilicé el algoritmo de línea de Xiaolin Wu para la rasterización. El programa genera el resultado en un archivo PGM llamado g.

Salida

#include<array>
#include<cmath>
#include<vector>
#include<string>
#include<fstream>
#include<algorithm>
#include<functional>
#define L for
#define A auto
#define E swap
#define F float
#define U using
U namespace std;
#define K vector
#define N <<"\n"
#define Z size_t
#define R return
#define B uint8_t
#define I uint32_t
#define P operator
#define W(V)<<V<<' '
#define Y template<Z C>
#define G(O)Y vc<C>P O(vc<C>v,F s){vc<C>o;L(Z i=0;i<C;++i){o\
[i]=v[i]O s;}R o;}Y vc<C>P O(vc<C>l, vc<C>r){vc<C>o;L(Z i=0;i<C;++i){o[i]=l[i]O r[i];}R o;}
Y U vc=array<F,C>;U v2=vc<2>;U v3=vc<3>;U v4=vc<4>;U m4=array<v4,4>;G(+)G(-)G(*)G(/)Y F d(
vc<C>a,vc<C>b){F o=0;L(Z i=0;i<C;++i){o+=a[i]*b[i];}R o;}Y vc<C>n(vc<C>v){R v/sqrt(d(v,v));
}v3 cr(v3 a,v3 b){R v3{a[1]*b[2]-b[1]*a[2],a[2]*b[0]-b[2]*a[0],a[0]*b[1]-b[0]*a[1]};}m4 P*(
m4 l,m4 r){R{l[0]*r[0][0]+l[1]*r[0][1]+l[2]*r[0][2]+l[3]*r[0][3],l[0]*r[1][0]+l[1]*r[1][1]+
l[2]*r[1][2]+l[3]*r[1][3],l[0]*r[2][0]+l[1]*r[2][1]+l[2]*r[2][2]+l[3]*r[2][3],l[0]*r[3][0]+
l[1]*r[3][1]+l[2]*r[3][2]+l[3]*r[3][3]};}v4 P*(m4 m,v4 v){R v4{m[0][0]*v[0]+m[1][0]*v[1]+m[
2][0]*v[2]+m[3][0]*v[3],m[0][1]*v[0]+m[1][1]*v[1]+m[2][1]*v[2]+m[3][1]*v[3],m[0][2]*v[0]+m[
1][2]*v[1]+m[2][2]*v[2]+m[3][2]*v[3],m[0][3]*v[0]+m[1][3]*v[1]+m[2][3]*v[2]+m[3][3]*v[3]};}
m4 at(v3 a,v3 b,v3 c){A f=n(b-a);A s=n(cr(f,c));A u=cr(s,f);A o=m4{1,0,0,0,0,1,0,0,0,0,1,0,
0,0,0,1};o[0][0]=s[0];o[1][0]=s[1];o[2][0]=s[2];o[0][1]=u[0];o[1][1]=u[1];o[2][1]=u[2];o[0]
[2]=-f[0];o[1][2]=-f[1];o[2][2]=-f[2];o[3][0]=-d(s,a);o[3][1]=-d(u,a);o[3][2]=d(f,a);R o;}
m4 pr(F f,F a,F b,F c){F t=tan(f*.5f);m4 o{};o[0][0]=1.f/(t*a);o[1][1]=1.f/t;o[2][3]=-1;o[2
][2]=c/(b-c);o[3][2]=-(c*b)/(c-b);R o;}F lr(F a,F b,F t){R fma(t,b,fma(-t,a,a));}F fp(F f){
R f<0?1-(f-floor(f)):f-floor(f);}F rf(F f){R 1-fp(f);}struct S{I w,h; K<F> f;S(I w,I h):w{w
},h{h},f(w*h){}F&P[](pair<I,I>c){static F z;z=0;Z i=c.first*w+c.second;R i<f.size()?f[i]:z;
}F*b(){R f.data();}Y vc<C>n(vc<C>v){v[0]=lr((F)w*.5f,(F)w,v[0]);v[1]=lr((F)h*.5f,(F)h,-v[1]
);R v;}};I xe(S&f,v2 v,bool s,F g,F c,F*q=0){I p=(I)round(v[0]);A ye=v[1]+g*(p-v[0]);A xd=
rf(v[0]+.5f);A x=p;A y=(I)ye;(s?f[{y,x}]:f[{x,y}])+=(rf(ye)*xd)*c;(s?f[{y+1,x}]:f[{x,y+1}])
+=(fp(ye)*xd)*c;if(q){*q=ye+g;}R x;}K<v4> g(F i,I r,function<v4(F,F)>f){K<v4>g;F p=i*.5f;F
q=1.f/r;L(Z zi=0;zi<r;++zi){F z=lr(-p,p,zi*q);L(Z h=0;h<r;++h){F x=lr(-p,p,h*q);g.push_back
(f(x,z));}}R g;}B xw(S&f,v2 b,v2 e,F c){E(b[0],b[1]);E(e[0],e[1]);A s=abs(e[1]-b[1])>abs
(e[0]-b[0]);if(s){E(b[0],b[1]);E(e[0],e[1]);}if(b[0]>e[0]){E(b[0],e[0]);E(b[1],e[1]);}F yi=
0;A d=e-b;A g=d[0]?d[1]/d[0]:1;A xB=xe(f,b,s,g,c,&yi);A xE=xe(f,e,s,g,c);L(I x=xB+1;x<xE;++
x){(s?f[{(I)yi,x}]:f[{x,(I)yi}])+=rf(yi)*c;(s?f[{(I)yi+1,x}]:f[{x,(I)yi+1}])+=fp(yi)*c;yi+=
g;}}v4 tp(S&s,m4 m,v4 v){v=m*v;R s.n(v/v[3]);}main(){F l=6;Z c=64;A J=g(l,c,[](F x,F z){R
v4{x,exp(-(pow(x,2)+pow(z,2))/(2*pow(0.75f,2))),z,1};});I w=1024;I h=w;S s(w,h);m4 m=pr(
1.0472f,(F)w/(F)h,3.5f,11.4f)*at({4.8f,3,4.8f},{0,0,0},{0,1,0});L(Z j=0;j<c;++j){L(Z i=0;i<
c;++i){Z id=j*c+i;A p=tp(s,m,J[id]);A dp=[&](Z o){A e=tp(s,m,J[id+o]);F v=(p[2]+e[2])*0.5f;
xw(s,{p[0],p[1]},{e[0],e[1]},1.f-v);};if(i<c-1){dp(1);}if(j<c-1){dp(c);}}}K<B> b(w*h);L(Z i
=0;i<b.size();++i){b[i]=(B)round((1-min(max(s.b()[i],0.f),1.f))*255);}ofstream f("g");f 
W("P2")N;f W(w)W(h)N;f W(255)N;L(I y=0;y<h;++y){L(I x=0;x<w;++x)f W((I)b[y*w+x]);f N;}R 0;}
  • l es la longitud de un lado de la cuadrícula en el espacio mundial.
  • c es el número de vértices a lo largo de cada borde de la cuadrícula.
  • La función que crea la cuadrícula se llama con una función que toma dos entradas, las coordenadas del espacio mundial xy z(+ y sube) del vértice, y devuelve la posición del espacio mundial del vértice.
  • w es el ancho de la pgm
  • h es la altura de la pgm
  • mes la vista / matriz de proyección. Los argumentos utilizados para crear mson ...
    • campo de visión en radianes
    • relación de aspecto de la pgm
    • cerca del plano de clip
    • plano de clip lejano
    • posición de la cámara
    • objetivo de la cámara
    • arriba vector

El procesador fácilmente podría tener más funciones, un mejor rendimiento y jugar mejor al golf, ¡pero me he divertido mucho!

Patrick Purcell
fuente
2
Wow, eso es increíble!
MD XF
1
Para nada ... ¡adelante!
Patrick Purcell
1
¡Ahí tienes, 133 bytes de descuento!
MD XF
1
¡Esto es fabuloso! Si pudieras decirme dónde aprendiste todo eso, ¡ sería genial !
HatsuPointerKun
1
@HatsuPointerKun Me alegro de que lo disfrutes! Este tutorial ... opengl-tutorial.org/beginners-tutorials/tutorial-3-matrices ... es un excelente lugar para comenzar.
Patrick Purcell el
9

Mathematica, 47 bytes

Plot3D[E^(-(x^2+y^2)/2/#^2),{x,-6,6},{y,-6,6}]&

toma como entrada σ

Entrada

[2]

salida
ingrese la descripción de la imagen aquí

-2 bytes gracias a LLlAMnYP

J42161217
fuente
1
Mathematica ganando? No hay sorpresas allí: P
MD XF
3
Guardar 2 bytes conE^(-(x^2+y^2)/2/#^2)
LLlAMnYP
6

R, 105 102 87 86 bytes

s=scan();plot3D::persp3D(z=sapply(x<-seq(-6,6,.1),function(y)exp(-(y^2+x^2)/(2*s^2))))

Toma Sigma de STDIN. Crea un vector de -6a 6en pasos de .1ambos xy y, luego crea una 121x121matriz tomando el producto externo de xy y. Esto es más corto que llamar matrixy especificar las dimensiones. La matriz ya está llena, pero está bien, porque estamos sobrescribiendo eso.

El forbucle-loop recorre los valores en x, haciendo uso de las operaciones vectorizadas en R, creando la matriz de densidad una fila a la vez.

(s)applyde nuevo es un método más corto para operaciones vectorizadas. Como el héroe que es, maneja la creación de la matriz por sí misma, ahorrando bastantes bytes.

ingrese la descripción de la imagen aquí

128 125 110 109 bytes, pero mucho más elegante:

Este diagrama es creado por el plotlypaquete. Lamentablemente, la especificación es un poco prolija, por lo que cuesta muchos bytes. Sin embargo, el resultado es realmente muy elegante. Recomiendo probarlo usted mismo.

s=scan();plotly::plot_ly(z=sapply(x<-seq(-6,6,.1),function(y)exp(-(y^2+x^2)/(2*s^2))),x=x,y=x,type="surface")

bla

JAD
fuente
Especifiqué en la pregunta que el gráfico no necesita tener números de línea, su segunda presentación está bien.
MD XF
Oh, debo haberme perdido eso. Cambié mis soluciones. Creo que la plotlytrama es lo suficientemente elegante como para garantizar que todavía se incluya aquí.
JAD
Bueno, ambos son mucho más elegantes que los míos : P
MD XF
Como solo lo usa suna vez, ¿podría hacer 2*scan()^2y quitar el s=scan();al principio? Ahorraría 3 bytes.
KSmarts
6

Applesoft BASIC, 930 783 782 727 719 702 695 637 bytes

-72 bytes y un programa de trabajo gracias a ceilingcat que detecta mi error y un algoritmo acortado

0TEXT:HOME:INPUTN:HGR:HCOLOR=3:W=279:H=159:L=W-100:Z=L/10:B=H-100:C=H-60:K=0.5:M=1/(2*3.14159265*N*N):FORI=0TO10STEPK:X=10*I+1:Y=10*I+B:HPLOTX,Y:FORJ=0TOL STEP1:O=10*J/L:D=ABS(5-I):E=ABS(5-O):R=(D*D+E*E)/(2*N*N):G=EXP(-R)*M:A=INT((C*G)/M):X=10*I+Z*O+1:Y=10*I+B-A:HPLOTTOX,Y:IF(I=0)GOTO4
1IF(J=L)GOTO3
2V=INT(J/10):IF((J/10)<>V)GOTO5
3D=ABS(5-I+K):E=ABS(5-O):R=(D*D+E*E)/(2*N*N):U=EXP(-R)/(2*3.14159*N*N):S=INT((C*U)/M):P=10*(I-K)+Z*O+1:Q=10*(I-K)+B-S:HPLOT TOP,Q:HPLOTX,Y
4IF(J=0)GOTO7:IF(I<10)GOTO5:IF(J=L)GOTO6:V=INT(J/10):IF((J/10)=V)GOTO6
5HCOLOR=0
6HPLOTTOX,10*I+B:HCOLOR=3:HPLOTX,Y
7NEXTJ:NEXTI:HPLOTW+1,H:HPLOTTO101,H:HPLOTTO0+1,H

Versión sin golf aquí.

Cuando se le da entrada 1:

entrada-1

Cuando se le da entrada 2:

entrada-2

MD XF
fuente
1
Esto muestra una vez más la superioridad de BASIC ....
Puede guardar algunos bytes más configurando una o más variables a un valor de uso frecuente, como 10. Además, sugiera reemplazar EXP(X)/(2*3.14159*S1*S1)conEXP(X)*M
ceilingcat el