genkuroki/MFmean by R.ipynb

## MFmean by R.ipynb
{
  "cells": [
    {
      "metadata": {
        "ExecuteTime": {
          "end_time": "2019-04-18T04:28:46.620993Z",
          "start_time": "2019-04-18T04:28:46.521Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "version",
      "execution_count": 1,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": "               _                           \nplatform       x86_64-w64-mingw32          \narch           x86_64                      \nos             mingw32                     \nsystem         x86_64, mingw32             \nstatus                                     \nmajor          3                           \nminor          5.3                         \nyear           2019                        \nmonth          03                          \nday            11                          \nsvn rev        76217                       \nlanguage       R                           \nversion.string R version 3.5.3 (2019-03-11)\nnickname       Great Truth                 "
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "library(ggplot2)\nlibrary(repr)\nlibrary(tidyr)\nlibrary(dqrng)",
      "execution_count": 2,
      "outputs": []
    },
    {
      "metadata": {
        "ExecuteTime": {
          "start_time": "2019-04-18T04:28:46.526Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "n <- 10^6\nN <- 1000\np <- 0.3\nt <- 2.0\n\nstart_time <- proc.time()\nMmean <- rep(0, n)\nFmean <- rep(0, n)\nfor (i in seq_len(n)) {\n  M <- rnorm(N)\n  F <- rnorm(N)\n  f <- F[F > t]\n  Mmean[i] <- mean(M[M > t])\n  Fmean[i] <- mean(f[1:(length(f) * p)])\n}\nproc.time() - start_time",
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": "   user  system elapsed \n 181.89    0.02  182.04 "
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {
        "ExecuteTime": {
          "start_time": "2019-04-18T04:28:46.527Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "options(repr.plot.width=5, repr.plot.height=3)\ndata.frame(Mmean, Fmean) %>%\n  tidyr::gather() %>%\n  ggplot(aes(x = value, fill = key)) +\n  geom_histogram(position = \"identity\", bins=50, alpha = 0.7) +\n  xlim(2, 3)",
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "stream",
          "text": "Warning message:\n\"Removed 532 rows containing non-finite values (stat_bin).\"Warning message:\n\"Removed 4 rows containing missing values (geom_bar).\"",
          "name": "stderr"
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": "plot without title",
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAFoCAMAAAC8KnXeAAAAUVBMVEUAAAAzMzNGzNBIztJJ\ntLVLtbdM0tZNTU1oaGh8fHyMjIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enr6+vw8PDy\n8vL0mZP2m5X6n5n///9X9oKHAAAACXBIWXMAABJ0AAASdAHeZh94AAAOv0lEQVR4nO2d7WKb\n1hJFT3Lb4g/Zsa2rpuX9H7RCAhshdAzDbEbAWj8kmezSMXsFEFZwKgEEpOgBYJ0gFkhALJCA\nWCABsUACYoEExAIJiAUSvMX6f5frJbch65F1btQIYq0u69yoEcRaXda5USOItbqsc6NGEGt1\nWedGjSDW6rLOjRpBrNVlnRs1gliryzo3agSxVpd1btQIYq0u69yoEcRaXda5USOItbqsc6NG\nEGt1WedGjaxcrJ8n/gqdYe6sc6NGEEs+w9xZ50aNIJZ8hrmzzo0aWatYP9sg1vwglm6GoKxz\no0YQSzdDUNa5USOIpZshKOvcqBHE0s0QlHVu1Ahi6WYIyjo3agSxdDMEZZ0bNYJYuhmCss6N\nGkEs3QxBWedGjSCWboagrHOjRhBLN0NQ1rlRI4ilmyEo69yoEcTSzRCUdW7UCGLpZgjKOjdq\nZBti1c8zzxCUdW7UCGLpZgjKOjdqBLF0MwRlnRs1gli6GYKyzo0aWeutIv/qI3qoLcEeSzdD\nUNa5USOIpZshKOvcqBHE0s0QlHVu1Ahi6WYIyjo3agSxdDMEZZ0bNYJYuhmCss6NGkEs3QxB\nWedGjSCWboagrHOjRhBLN0NQ1rlRI6sT62cPiDU/iCWYITbr3KgRxBLMEJt1btQIYglmiM06\nN2oEsQQzxGadGzWCWIIZYrPOjRpBLMEMsVnnRo0glmCG2Kxzo0YQSzBDbNa5USOIJZghNuvc\nqBHEEswQm3Vu1AhiCWaIzTo3agSxBDPEZp0bNYJYghlis86NGkEswQyxWedGjSCWYIbYrHOj\nRhBLMENs1rlRI4glmCE269yokU2J9XOAYPcpy5isc6NGEEswQ2zWuVEjiCWYITbr3KgRxBLM\nEJt1btTI92IVR4Y8nwncoGcQS2CJgW/FKuqH755rAjfoGcQSWGIAsQQzxGYFlhhALMEMsVmB\nJQYcxfpRoZpzML03H+VepLMzSKyCPdaCsgpNxjPoXSGHwiVlFZqMZ9h1LMRaUNZbERucvAtm\niM0KLDGAWIIZYrMCSwxw5V0wQ2xWYIkBflYomCE269yoEcQSzBCbdW7UCGIJZojNOjdqBLEE\nM8RmnRs1gliCGWKzzo0aQSzBDLFZ50aNIJZghtisc6NGEEswQ2zWuVEjiCWYITbr3KgRxBLM\nEJt1btQIYglmiM06N2oEsQQzxGadGzWCWIIZYrPOjRpBLMEMsVnnRo2sR6w+oxArDMSKmleW\ndW7UCGJFzSvLOjdqBLGi5pVlnRs1glhR88qyzo0aQayoeWVZ50aNIFbUvLKsc6NGECtqXlnW\nuVEjiBU1ryzr3KgRxIqaV5Z1btQIYkXNK8s6N2oEsaLmlWWdGzWySbFyft2nLGOyzo0aQayo\neWVZ50aNeIsVx5B7RHLHyNlgjxU1ryzr3KgRxIqaV5Z1btQIYkXNK8s6N2oEsaLmlWWzdc92\nTo1YUfPKstm6EWt0FrGaJbm6EWt0FrGaJbm6q7536a0sD88pPR/KQ3qoltdPjiBW1LyybLbu\nVHm1O74q0pGH2rLyV3pxFgGxouaVZbN1p6NJlUMvlV279Fp+pMfjl0/p3VkExIqaV5bN1p3O\ne6jy4VR8eqqc+ji+KHL/lQXEippXls3WfeS1eVFRHndZT+Vbenb2ALHC5pVls3WnlyLty7ZY\nx53Xvt6NeYJYUfPKstm60/E8/Xj8qw+FJ97SrvC/CoFYUfPKstm6j30/pl/1W8NfpxP3o2Tn\nZ1cQK2peWTZb9+mkqjiUh9Plhuq8vdplVao5g1hR88qy2bqrvl+qU/X983E/9d4s3TtrgFhx\n88qyozt7d7/sXiJW3Lyy7OjOHv3fEyJW3Lyy7FgDFKfuiBU3ryw7srHidPnBHcSKmleWdW7U\nSEes5vM6hfVnR3EbFLGaJXdBW6witTCuL26DIlaz5C5oC/Ta8urVuL64DYpYzZK74Mah0Ezc\nBkWsZsldwMl71LyyrHOjRrpi7QrOseaZV5adqoQPHYF2nLzPNa8sO9kJFzoCFeaz9pq4DYpY\nzZIM/9zm72nFd+HkPWpeWTZXT5hYT+kwbX1xG7S25Y8ziNVLmFj74rH70ZyiOF+G/+75TNwG\nvRTrwi/Eagg8FHZP3ov64bvnmrgNiljNkgyINT77Rx+IdUmYWNcU7WfEWkA21+ZdiXU+h/pe\nrB8VvtON4c8M3IO05p4OhbVVy9xjnfmTPVbNPYlVP9ytWI0miNUsyTBMrGk/fOkV68z+8eum\nNojlPK88m6t7oFg2ly7oX8fh63ZJiOU8rzybqztarHJBh8KMUN+I1evXfcoyJpure6RYqTwd\nEevDYrp8TmXucNn/J79a90u69yvviNVd4ihW+nyol7af001/rv/g89x9lxsvw/wbdIhY9TNi\njTl5L78cuhCq7Pv6in6xCqtXiBWfzdUz+lBYtvdQl8fCUWJNZv4NiljdJRkmidVeiliI1cZD\nrN5D4xXdPzjsHlJ62Jk/lTX/BkWs7pIMk8Wqn0cfCvf1v6UorDdMmn+DIlZ3SYYpYrUvN4ze\nYz2n6oN++0fzXXTn36CI1V2SYZhYHtz4zPuC/pUOYnWXZECsAdnfZxCruyRDmFhLOhSOF6sB\nsWYXa0kn74jVn83VEybWki43IFZ/NldPnFhTmXGDIlZ/NlfP3xl8RUAsxJpFrKfzZa8HzrHc\n550rO00IL67vNnNayrtC/3nnyubq+d9ttHus4vybNj/WeR0LscLEWvcFUsQKE+up+s3m5WFn\n/mUFM25QxOrP5uoJE+vzAumHcX0zblC7WJd+zTXvXNlcPWFiNRdIzb9mbMYNilj92Vw9cWJN\nZY4NevLp39+I1Z/N1YNYORArm83Vg1g5ECubzdWDWDkQK5vN1TNMrNR5toBYiHUtVmo/2UAs\nxEKsGsTKZnP1DD0UpubRfFsQb7Hm4N8LcreIzLPBO0aOF8t6WxD2WOyx8nus8lqosu/rDlsW\nq0Y+78zZXD1Dxfq0ynxbEMRCrOFita1BLMT6wkesIbcFQSzE6hPreqc08rYgiIVYg8Uac1sQ\nxOq8ObxPWcZkc/UME8sDxEIsxKpBrGw2Vw9i9fK7BWLdyubqQaxeEGtINlcPYvWCWEOyuXoQ\nqxfEGpLN1cNNQXpBrCFZ50aNLFisAUYNEqtGMG9M1rlRI4iFWBIQC7EkIBZiSUAsxJKAWIgl\nAbEQSwJiIZYExEIsCYiFWBIQC7EkLEOs39cg1q2sc6NGEKuhXu99yjIm69yoEcRCLAmIhVgS\nEAuxJCAWYkn4XqziyJDnM5oNilgjsgJLDHwrVlE/fPdco9mgiDUiK7DEAGIhloRh51iINW7e\n0Ky3IjYcxfpRIRmyc9fRyTcf7aVeq+Yb2B4DxDqfnG9lj9X8s7JJ84Zm/SWxwKGwA2L5gFiI\nJYF3hYglAbEQSwJX3hFLAj8rRCwJiIVYEu5brB6hEOu7rHOjRhDrhlhD/EKs2yAWYklALMSS\ngFiIJQGxEEsCYiGWBMRCLAmIhVgSlijWaFkQa34QC7EkIFZerKxfiHUbxEIsCYiFWBIQC7Ek\nIBZiSbhTsfqNQqwhWedGjSDWjSxiTQOxEEsCYg0Sq98vxLoNYiGWBMRCLAmIhVgSEAuxJCAW\nYklALMSS4C2WE733hhTdIzJ758heorfOEmCPNWaP1dlxsce6DWIhlgTEQiwJ9yZW1ijEGpJ1\nbtQIYiGWBMRCLAmIlc9m/UKs2yAWYklALMSSsCixnGRBrBlALMSSgFh2sZrfFGYUQJV1btQI\nYiGWBMRCLAmIhVgS7kasVlFLE+v3EL8QaxrmjYRYXlnnRo0g1qAsYo0FsRBLAmK5iZX3C7Gm\nYd5IiOWVdW7UCGIhlgTE8hbrhl+INY3xG2lIURJZEEsJYo3JItZgEEskVtcvxJrG+I2EWOO3\nWTbr3KgRxNKK9SkYYk1j/Ebahlh11irLmKxzo0bixLq58e9YrBrE+p7pYhVHvr4avpGGiDWj\nLGOydrFqEOt7is+HE5sQq8Yu1iWI1cd4sUZs/PlliRDryi/EKlVihckyITtVrP4sYpU/Kqau\nDlZC3DnWoL99ZMdnpzbqA2KtLju1UR8Qa3XZqY36gFiry05t1AfEWl12aqM+xF15d9+gZOsl\nd0H8D6HJOmedGzWCWKvLOjdqBLFWl3Vu1AhirS7r3KgRxFpd1rlRI4i1uqxzo0YQa3VZ50aN\nyH8RpurzDqxXu96pIBbrlYBYrFcCYrFeCXf6y8Zh6SAWSEAskIBYIAGxQAJigQSJWO0PlTav\nLz9o6rveqSvurrco3ect2uv1nre77B5QiNX+GHzzuvPReLf1Tl1nd71X6/dab/PaeV7P7esK\nYs0nllf5mxXrhOob7+wB3LZlcfEkmdfp78HFehGrtcxjvc1rj1OW9no/T4W6/6+J6/1cpfe8\n9cNmxOr+zXc9DLRee69XOG9x+TBprYVqXjeWLtbVQof1KsW6/p9NW/HWxLrs3O8b7/XKeb2i\neftfTl3zxsTqnli4HQKuXnuvdwnzbvddYc+Jhcs33rO/cnmn1fPOzX3ellguq92iWEX7rZXj\nleH2ei+uZC9g3i8THK/oO8/rCj8rBAmIBRIQCyQgFkhALJCAWCABsUACYoEExAIJiHVBYns4\nwYa8ALG8YENegFhebHdDHtLD6fkhfZTvTykVu/Is1lmu0+PhOaXnQ+CQy2W7YpVPaX983B/9\neksndldiFdXih9gxF8qGxXqrTCp36e240/pVlh+NVF9ivVSJXXoNHnSRbFis8uH8eabq5f7t\n5fFarIfzq6fAGRfLlsV6Te/le3o5vno8Hwu7YqXULIexbHmjHdLz8UB3PDd/Tg+vb3vE8mTT\nG+057U/HufM7wAux9l+HQrCw6U33ftwbvZeVSu/l4fMcqzieyp+/2lUn77/SY/ScS2TTYh13\nSadrCbvUPsc6ffVSvTqcLjekj+gxl8i2xXqtrjOU1TExPb5/HgZ3xfGM/nxAPP1B6IhLZdti\ngQzEAgmIBRIQCyQgFkhALJCAWCABsUACYoEExAIJiAUSEAsk/Af0ls+gU/JpKgAAAABJRU5E\nrkJggg=="
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {
        "ExecuteTime": {
          "start_time": "2019-04-18T04:28:46.526Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "n <- 10^6\nN <- 1000\np <- 0.3\nt <- 2.0\n\nstart_time <- proc.time()\nMmean <- rep(0, n)\nFmean <- rep(0, n)\nfor (i in seq_len(n)) {\n  M <- dqrnorm(N)\n  F <- dqrnorm(N)\n  f <- F[F > t]\n  Mmean[i] <- mean(M[M > t])\n  Fmean[i] <- mean(f[1:(length(f) * p)])\n}\nproc.time() - start_time",
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": "   user  system elapsed \n  37.12    0.09   37.25 "
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {
        "ExecuteTime": {
          "start_time": "2019-04-18T04:28:46.527Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "options(repr.plot.width=5, repr.plot.height=3)\ndata.frame(Mmean, Fmean) %>%\n  tidyr::gather() %>%\n  ggplot(aes(x = value, fill = key)) +\n  geom_histogram(position = \"identity\", bins=50, alpha = 0.7) +\n  xlim(2, 3)",
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "stream",
          "text": "Warning message:\n\"Removed 538 rows containing non-finite values (stat_bin).\"Warning message:\n\"Removed 4 rows containing missing values (geom_bar).\"",
          "name": "stderr"
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": "plot without title",
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAFoCAMAAAC8KnXeAAAAUVBMVEUAAAAzMzNGzNBIztJJ\ntLVLtbdM0tZNTU1oaGh8fHyMjIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enr6+vw8PDy\n8vL0mZP2m5X6n5n///9X9oKHAAAACXBIWXMAABJ0AAASdAHeZh94AAAOHUlEQVR4nO2dbWOi\nuhZGM3PPOdjWvno9vcP//6FXFCsKSRPybCLNWh/E0mcy282aEKlDXQtggCtdAPxMEAtMQCww\nAbHAhBixmtPjgZgtQBslVu9N//DdFqDje7GaFrEgmfhTYf8MsSAGA7F+dajrhJWRJFbTRs9Y\n/71lvMcL0axhhXZkgFj3GK1MrKhFfI+0VdVF6xKruTwilm20KrGGeiGWbbQmsZqmv7Qed+Vd\n2qrqopWINQdpq6qLIpYXaauqiyKWF2mrqosilhdpq6qLIpYXaauqiyKWF2mrqosilhdpqzKj\nv4/8U66A5ChieZG2KjOKWAVArAUKSI4ilhdpqzKjiFUAxFqggOQoYnmRtiozilgF+NFi/R6C\nWIuCWJYFzI0ilhdpq+ZGEasciGVZwNwoYnmRtmpuFLHKgViWBcyNIpYXaavmRhGrHIhlWcDc\nKGJ5kbZqbhSxyoFYlgXMjSKWF2mr5kYRqxwVidVvlyxgbhSxvEhbNTeKWOVALMsC5kYRy4u0\nVXOjiFUOxLIsYG4UsbxIWzU3iljl+NG/QOCfKUoXVQnMWJYFzI0yY3mRtmpuFLHKgViWBcyN\nIpYXaavmRhGrHIhlWcDcKGJ5kbZqbhSxyoFYlgXMjSKWF2mr5kYRqxyIZVnA3ChieZG2am4U\nscrxM8X6PQFiLQpimRSQGUUsL9JWJUcRqziIZVJAZhSxvEhblRxFrOIglkkBmVHE8iJtVXIU\nsYqDWCYFZEYRy4u0VclRxCoOYpkUkBlFLC/SViVHEas4iGVSQGYUsbxIW5UcRaziIJZJAZlR\nxPIibVVyFLGKU59YEX4hVj6IZVJAZhSxvEhblRxFrOIglkkBmVHE8iJtVXIUsYoTL1ZzIGZ7\nQtqq5ChiFSdarKZ/+G7bI21VchSxioNYJgVkRhELsUyiiDUl1q8OdZ1pTN4jkjtGLkmSWA0z\n1jLRmmas7l0fp8KFolWJ1YFYy0RrEovF+4JRxEIsk2hNYnHlfcFoVWKlIW1VchSxioNYJgVk\nRhHLi7RVyVHEKg5imRSQGUUsL9JWJUcRqziIZVJAZhSxvEhblRxFrOIglkkBmVHE8iJtVXIU\nsYqDWCYFZEYRy4u0VclRxCoOYpkUkBlFLC/SViVHEas4iGVSQGYUsbxIWxUfnTIKsUqAWNoC\nNFHE8iJtVXwUsRDLpKuIhVgmXUUsxDLpKmIhlklXEQuxTLqKWIhl0lXEQiyTriIWYpl0FbEQ\ny6SriIVYJl1FLMQy6Wq8WCG/ECsfK7EKEXGPSO4YuQjMWNoCNFFmLC/SVsVHEQuxTLqKWGGx\n3GJLH8TSFqCJIpYXaavio4iFWCZdRawIsbburW33T8497du923T7+40QxNIWoImairV128Oz\nxh3Y9Ja1r+5ZLABiaQvQRC3F2h4deu7s2rqX9sM9HL58dO9iARBLW4AmaijWaYZqN8cD7x47\npz4OT5rQn5oDYmkL0EQNxXKHWap/0tEepqzH9s09qQVALG0BmqihWM+N27VDsQ6T166fxpQg\nlrYATdRyjfXanf/6U+GRN7dt9BoglrYATdT0XeGDe+3fGr4eF+4HyU5bKYilLUATNRXrwzX7\ndn+83NCt27spq1NNDGJpC9BEbS+QPndL9d3TYZ56P+/dyQVALG0Bmuiyn254l192bxFLXYAm\nuqxYD/r3hIilLkATXVIsZ7F0Ryx1AZrokmI1x8sPchBLW4Amuuyp0IQbsc6f12lyf3YkbVV8\nFLHuUazGDcgcV9qq+Chi3aNYLwOvXjLHlbYqPopY9yhWK/zoqrRV8VHEulOxZEhbFR9FrHsV\na9uwxkIsATcCbVm8L1ar1bCZR07EjUBN9qq9R9qq+ChihcX6n59/NQf+DIt3bQGa6M8T69Ht\nNeNKWxUf7XX56wRi3VJMrF3z4PtoTtOcLsd/tz0hbVV89FqsK78Qqy16KvQt3pv+4bttj7RV\n8VHEQiyTrv41BWJ9UUwsP81wi1i20brEOq2hvhfrV4e2ymj+DsCtIu/zVNhbtcIZ68TfzFj3\nKVb/cKdiTa7aEeuWOLEUP3yZ/NO7h/FNbRBrwWhpsdIcmmR6jP34dkmItWD0x4o18aOdOxcr\nYBRiXUgUy7XHM2J/WnTXW9eGTpfT33mduF/SfV95nyvWpF+IdX7ivh76vcOt8/oz/sbX2n0b\nKi8Caasiooh12REgZfHeXhy6Eqqd+nrEtFhNrlf3KFa/RayUU2E7nKGuz4VJYsmQtioiiliX\nHQGyxBruRSzEGqIQa/LUOOL2G/vtxrnNNvtTWdJWRUQR67IjQLZY/Tb5VLjr/y9Fk3vDJGmr\nIqKIddkRIEes4eWG5BnryXUf9Ns9ZN9FV9qqiChiXXYEiBNLgecz76v7XzqIddkRALFSo4h1\n2RGgmFjrOxV+nkCsy44AxcRa3+I9WaxJv5apdZlhQ4elmFjru9yAWKMdAcqJpULaqlAUsUY7\nAvwbQCsAYiHWImI9ni57bVhjGdS6zLCZR07E+G4zx70/+F1h1WL9x4/tjNWcftPmx3quYyHW\naEeAYmKt7wIpYo12BCgm1mP3m83b/Tb7lxVIWxWKItZoR4BiYn1dIP3IHFfaqlAUsUY7AhQT\n63yBNPvXjElbFYrOFevar2VqXWbY0GEpJ5YKaasmOfr05xOxRjsCINb3IJYnGjosiPU9iOWJ\nhg4LYn0PYnmiocMSJ5a72c4BsRBrLJYbbuaBWIiFWFcglicaOiyxp0J3fpx9WxDEQqygWHNv\nC2Illj1/rgjdfDTixqSlX8xSzJix2nm3BWHGYsaaEOvLqtm3BUEsxIoXaxBgxgqL1WNca2r0\nzsWKuS0IYiHWlFjjSSnxtiCIhVjRYqXcFgSxEGsklgLEOnK9hkesfFYo1ucAxEIsWVcRKxwN\nHRbECoBY4WjosCBWAMQKR0OHBbECIFY4Gjos3BQkAGKFo0ZHNBHEGiKvdWYUsbxIW3XNtVgR\nuiBWARALsUxALMQyAbEQywTEQiwTEAuxTEAsxDIBsRDLBMRCLBMQC7FMQCzEMgGxhvSjIlY+\naxLrcwxiIVZ+VxErKmp0RBNBLMQyAbEQy4R4sZoDMdsT0lb1IFZUVGhHBtFiNf3Dd9seaat6\nECsqKrQjA8RCLBPS1liIZfCy1MOq1MjDQKxfHdIie/5MMfcekZP0Y1oUXxsJYp0W5zXMWP3/\n/cmpNS/KjIVY+S9LPaxKjTwQawLEyod3hYhlAmIhlglceUcsE/hZIWKZsAqxJoxCLG/U6Igm\ngliIZQJiIZYJiOUX6zPCL8TygFiIZQJiIZYJiIVYJiAWYpmAWIhlAmIhlgmIhVgmrFesVFsQ\na1EQC7FMQKxvxQr6hVgeEAuxTEAsxDIBsRDLhPsWa9IoxApHjY5oIojljyJWBoiFWCYgFmKZ\ngFixYk37hVgeEAuxTEAsxDIBsRDLBMRCLBMQC7FMsBJLw+TNIS3uERm6ceQ1pVuyFpixEmes\nm5mLGcsDYiGWCYiFWCYgFmKZcKdiBY8pYgWjRkc0EcRCLBMQ69tosArE8oBYiGXCCsWS2IJY\nxiAWYpmAWIhlAmLliWVzH1zE8pLZKsSaHzU6ookgFmKZgFiIZcK9iTU4WIg1L2p0RBNBLIlY\nnxF+IZaCua26Q7GCfiGWB8RCLBMQSylW0C/EUjC3VYiVmkWsKBArNYtYUfwIsab9QiwFya2K\nOFgGtiCWFYiVGEWsOBDLTqwbvxBLQXKrECulW6Go0RFNBLHMxTr7hVgKklu1FrHOJIvVR2fa\nkhA1OqKJlBfLdwQQC7GONAcuX0nFWsyWhOhcsXoQK5bm6+HIjxerZ65Y1yCWn3Sx4o/A0rYs\nLtbIL8T6wkSsQrYs8ZnAMIMsYp0ffnWohoWVUn6NFfFvsLoop8ILiCWMItYFxBJGEesCYgmj\niHUBsYRRxBow98q7tqs/I4pYXqStqi6KWF6kraouilhepK2qLopYXqStqi6KWF6kraouilhe\npK2qLopYXqStqi6KWPHYfN6BUa2GzQWxVj4qYjGqyaiIxagmo1YvFtQFYoEJiAUmIBaYgFhg\nAmKBCaZiDT9Uen5+/UFT5ah5w96O2rTiWpvhqEYtyBtViqVYw4/Bn5/ffDReNGreiLejjkbX\njHp+nn/4bRqrBbHGo9qJJTr6tYt1xOb138wCooY2VxuDWjUT1tWwiNVaiqVZYJxrPS+Gbv+m\nrFG/BhQthioX6/Zfv+b1Ty5bZKOa1dpcP+QM2tgUq+NniDXamT2qnVjjvypr3ErFuj7qqtc/\n6ZV0VJNap59mDlynWLfLC83rH3ulHfXua63+XeHE8kLV1uvngndbE+/dxLUOxFLUWrFYzfDt\nlewC8XDUq6vZd17rRYX8d4VfBbbCYrXws0IwAbHABMQCExALTEAsMAGxwATEAhMQC0xALDAB\nsUY4eiKAJo5ALAU0cQRiKai7iXu3OW437qN9f3Su2bYnsU5yHR/3T8497QsWuU7qFqt9dLvD\n4+7g15s7sh2J1XS7N2XLXCGVi/XWmdRu3dth0npt24+zVBexnrvE1r0ULnR1VC5Wuzl9pKl7\nunt7fhiLtTk9eyxY4yqpXawX996+u+fDs4fTufBWLOfO+yGF2hu2d0+HE91hbf7kNi9vO8RS\nUX3DntzueJ47vQO8Emt3ORVCKtW37f0wG723nUrv7f5rjdUclvKnr7bd4v3VPZSuc21UL9Zh\nSjpeS9i64Rrr+NVz92x/vNzgPkqXuTYQ66W7ztB250T38P51Gtw2hxX96YR4/EbREtcIYoEJ\niAUmIBaYgFhgAmKBCYgFJiAWmIBYYAJigQmIBSYgFpiAWGDC/wGwkyHZzfyOSgAAAABJRU5E\nrkJggg=="
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "simpi <- function (n) {\n    s <- 0\n    for (i in seq_len(n)) {\n        x <- runif(2)\n        s <- s + ifelse(x[1]^2 + x[2]^2 < 1.0, 1, 0)\n    }\n    4*s/n\n}\n\nstart_time <- proc.time()\nsimpi(10^7)\nproc.time() - start_time",
      "execution_count": 17,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/html": "3.1403876",
            "text/markdown": "3.1403876",
            "text/latex": "3.1403876",
            "text/plain": "[1] 3.140388"
          },
          "metadata": {}
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": "   user  system elapsed \n  31.83    0.01   31.85 "
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "simpi_dqr <- function (n) {\n    s <- 0\n    for (i in seq_len(n)) {\n        x <- dqrunif(2)\n        s <- s + ifelse(x[1]^2 + x[2]^2 < 1.0, 1, 0)\n    }\n    4*s/n\n}\n\nstart_time <- proc.time()\nsimpi_dqr(10^7)\nproc.time() - start_time",
      "execution_count": 18,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/html": "3.1424212",
            "text/markdown": "3.1424212",
            "text/latex": "3.1424212",
            "text/plain": "[1] 3.142421"
          },
          "metadata": {}
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": "   user  system elapsed \n  22.68    0.00   22.72 "
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "simpi_dqr1000 <- function (n) {\n    s <- 0\n    N = 1000\n    m = n/N\n    for (i in seq_len(m)) {\n        x <- dqrunif(N)\n        y <- dqrunif(N)\n        for (j in seq_len(N)) {\n            s <- s + ifelse(x[j]^2 + y[j]^2 < 1.0, 1, 0)\n        }\n    }\n    4*s/n\n}\n\nstart_time <- proc.time()\nsimpi_dqr1000(10^7)\nproc.time() - start_time",
      "execution_count": 25,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/html": "3.1411412",
            "text/markdown": "3.1411412",
            "text/latex": "3.1411412",
            "text/plain": "[1] 3.141141"
          },
          "metadata": {}
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": "   user  system elapsed \n  15.62    0.00   15.62 "
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "simpi_dqr_sum <- function (n) {\n    x <- dqrunif(n)\n    y <- dqrunif(n)\n    4*sum(x^2 + y^2 < 1.0)/n\n}\n\nstart_time <- proc.time()\nsimpi_dqr_sum(10^7)\nproc.time() - start_time",
      "execution_count": 39,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/html": "3.1406348",
            "text/markdown": "3.1406348",
            "text/latex": "3.1406348",
            "text/plain": "[1] 3.140635"
          },
          "metadata": {}
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": "   user  system elapsed \n   0.19    0.16    0.41 "
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "start_time <- proc.time()\nsimpi_dqr_mean(10^8)\nproc.time() - start_time",
      "execution_count": 40,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/html": "3.14144876",
            "text/markdown": "3.14144876",
            "text/latex": "3.14144876",
            "text/plain": "[1] 3.141449"
          },
          "metadata": {}
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": "   user  system elapsed \n   1.95    2.06    4.02 "
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "start_time <- proc.time()\nsimpi_dqr_mean(10^9)\nproc.time() - start_time",
      "execution_count": 41,
      "outputs": [
        {
          "output_type": "error",
          "ename": "ERROR",
          "evalue": "Error:  サイズ 7.5 Gb のベクトルを割り当てることができません \n",
          "traceback": [
            "Error:  サイズ 7.5 Gb のベクトルを割り当てることができません \nTraceback:\n",
            "1. simpi_dqr_mean(10^9)",
            "2. mean(x^2 + y^2 < 1)   # at line 4 of file <text>"
          ]
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "",
      "execution_count": null,
      "outputs": []
    }
  ],
  "metadata": {
    "kernelspec": {
      "name": "ir",
      "display_name": "R",
      "language": "R"
    },
    "language_info": {
      "name": "R",
      "codemirror_mode": "r",
      "pygments_lexer": "r",
      "mimetype": "text/x-r-source",
      "file_extension": ".r",
      "version": "3.5.3"
    },
    "toc": {
      "nav_menu": {},
      "number_sections": true,
      "sideBar": true,
      "skip_h1_title": false,
      "title_cell": "Table of Contents",
      "title_sidebar": "Contents",
      "toc_cell": false,
      "toc_position": {},
      "toc_section_display": true,
      "toc_window_display": false
    },
    "gist": {
      "id": "9b2f093157d8ee9e4587fbf06eb0ba5e",
      "data": {
        "description": "MFmean by R.ipynb",
        "public": true
      }
    },
    "_draft": {
      "nbviewer_url": "https://gist.github.com/9b2f093157d8ee9e4587fbf06eb0ba5e"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 2
}
	{
	"cells": [
	{
	"metadata": {
	"ExecuteTime": {
	"end_time": "2019-04-18T04:28:46.620993Z",
	"start_time": "2019-04-18T04:28:46.521Z"
	},
	"trusted": true
	},
	"cell_type": "code",
	"source": "version",
	"execution_count": 1,
	"outputs": [
	{
	"output_type": "display_data",
	"data": {
	"text/plain": " _ \nplatform x86_64-w64-mingw32 \narch x86_64 \nos mingw32 \nsystem x86_64, mingw32 \nstatus \nmajor 3 \nminor 5.3 \nyear 2019 \nmonth 03 \nday 11 \nsvn rev 76217 \nlanguage R \nversion.string R version 3.5.3 (2019-03-11)\nnickname Great Truth "
	},
	"metadata": {}
	}
	]
	},
	{
	"metadata": {
	"trusted": true
	},
	"cell_type": "code",
	"source": "library(ggplot2)\nlibrary(repr)\nlibrary(tidyr)\nlibrary(dqrng)",
	"execution_count": 2,
	"outputs": []
	},
	{
	"metadata": {
	"ExecuteTime": {
	"start_time": "2019-04-18T04:28:46.526Z"
	},
	"trusted": true
	},
	"cell_type": "code",
	"source": "n <- 10^6\nN <- 1000\np <- 0.3\nt <- 2.0\n\nstart_time <- proc.time()\nMmean <- rep(0, n)\nFmean <- rep(0, n)\nfor (i in seq_len(n)) {\n M <- rnorm(N)\n F <- rnorm(N)\n f <- F[F > t]\n Mmean[i] <- mean(M[M > t])\n Fmean[i] <- mean(f[1:(length(f) * p)])\n}\nproc.time() - start_time",
	"execution_count": 3,
	"outputs": [
	{
	"output_type": "display_data",
	"data": {
	"text/plain": " user system elapsed \n 181.89 0.02 182.04 "
	},
	"metadata": {}
	}
	]
	},
	{
	"metadata": {
	"ExecuteTime": {
	"start_time": "2019-04-18T04:28:46.527Z"
	},
	"trusted": true
	},
	"cell_type": "code",
	"source": "options(repr.plot.width=5, repr.plot.height=3)\ndata.frame(Mmean, Fmean) %>%\n tidyr::gather() %>%\n ggplot(aes(x = value, fill = key)) +\n geom_histogram(position = \"identity\", bins=50, alpha = 0.7) +\n xlim(2, 3)",
	"execution_count": 4,
	"outputs": [
	{
	"output_type": "stream",
	"text": "Warning message:\n\"Removed 532 rows containing non-finite values (stat_bin).\"Warning message:\n\"Removed 4 rows containing missing values (geom_bar).\"",
	"name": "stderr"
	},
	{
	"output_type": "display_data",
	"data": {
	"text/plain": "plot without title",
	"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAFoCAMAAAC8KnXeAAAAUVBMVEUAAAAzMzNGzNBIztJJ\ntLVLtbdM0tZNTU1oaGh8fHyMjIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enr6+vw8PDy\n8vL0mZP2m5X6n5n///9X9oKHAAAACXBIWXMAABJ0AAASdAHeZh94AAAOv0lEQVR4nO2d7WKb\n1hJFT3Lb4g/Zsa2rpuX9H7RCAhshdAzDbEbAWj8kmezSMXsFEFZwKgEEpOgBYJ0gFkhALJCA\nWCABsUACYoEExAIJiAUSvMX6f5frJbch65F1btQIYq0u69yoEcRaXda5USOItbqsc6NGEGt1\nWedGjSDW6rLOjRpBrNVlnRs1gliryzo3agSxVpd1btQIYq0u69yoEcRaXda5USOItbqsc6NG\nEGt1WedGjaxcrJ8n/gqdYe6sc6NGEEs+w9xZ50aNIJZ8hrmzzo0aWatYP9sg1vwglm6GoKxz\no0YQSzdDUNa5USOIpZshKOvcqBHE0s0QlHVu1Ahi6WYIyjo3agSxdDMEZZ0bNYJYuhmCss6N\nGkEs3QxBWedGjSCWboagrHOjRhBLN0NQ1rlRI4ilmyEo69yoEcTSzRCUdW7UCGLpZgjKOjdq\nZBti1c8zzxCUdW7UCGLpZgjKOjdqBLF0MwRlnRs1gli6GYKyzo0aWeutIv/qI3qoLcEeSzdD\nUNa5USOIpZshKOvcqBHE0s0QlHVu1Ahi6WYIyjo3agSxdDMEZZ0bNYJYuhmCss6NGkEs3QxB\nWedGjSCWboagrHOjRhBLN0NQ1rlRI6sT62cPiDU/iCWYITbr3KgRxBLMEJt1btQIYglmiM06\nN2oEsQQzxGadGzWCWIIZYrPOjRpBLMEMsVnnRo0glmCG2Kxzo0YQSzBDbNa5USOIJZghNuvc\nqBHEEswQm3Vu1AhiCWaIzTo3agSxBDPEZp0bNYJYghlis86NGkEswQyxWedGjSCWYIbYrHOj\nRhBLMENs1rlRI4glmCE269yokU2J9XOAYPcpy5isc6NGEEswQ2zWuVEjiCWYITbr3KgRxBLM\nEJt1btTI92IVR4Y8nwncoGcQS2CJgW/FKuqH755rAjfoGcQSWGIAsQQzxGYFlhhALMEMsVmB\nJQYcxfpRoZpzML03H+VepLMzSKyCPdaCsgpNxjPoXSGHwiVlFZqMZ9h1LMRaUNZbERucvAtm\niM0KLDGAWIIZYrMCSwxw5V0wQ2xWYIkBflYomCE269yoEcQSzBCbdW7UCGIJZojNOjdqBLEE\nM8RmnRs1gliCGWKzzo0aQSzBDLFZ50aNIJZghtisc6NGEEswQ2zWuVEjiCWYITbr3KgRxBLM\nEJt1btQIYglmiM06N2oEsQQzxGadGzWCWIIZYrPOjRpBLMEMsVnnRo2sR6w+oxArDMSKmleW\ndW7UCGJFzSvLOjdqBLGi5pVlnRs1glhR88qyzo0aQayoeWVZ50aNIFbUvLKsc6NGECtqXlnW\nuVEjiBU1ryzr3KgRxIqaV5Z1btQIYkXNK8s6N2oEsaLmlWWdGzWySbFyft2nLGOyzo0aQayo\neWVZ50aNeIsVx5B7RHLHyNlgjxU1ryzr3KgRxIqaV5Z1btQIYkXNK8s6N2oEsaLmlWWzdc92\nTo1YUfPKstm6EWt0FrGaJbm6EWt0FrGaJbm6q7536a0sD88pPR/KQ3qoltdPjiBW1LyybLbu\nVHm1O74q0pGH2rLyV3pxFgGxouaVZbN1p6NJlUMvlV279Fp+pMfjl0/p3VkExIqaV5bN1p3O\ne6jy4VR8eqqc+ji+KHL/lQXEippXls3WfeS1eVFRHndZT+Vbenb2ALHC5pVls3WnlyLty7ZY\nx53Xvt6NeYJYUfPKstm60/E8/Xj8qw+FJ97SrvC/CoFYUfPKstm6j30/pl/1W8NfpxP3o2Tn\nZ1cQK2peWTZb9+mkqjiUh9Plhuq8vdplVao5g1hR88qy2bqrvl+qU/X983E/9d4s3TtrgFhx\n88qyozt7d7/sXiJW3Lyy7OjOHv3fEyJW3Lyy7FgDFKfuiBU3ryw7srHidPnBHcSKmleWdW7U\nSEes5vM6hfVnR3EbFLGaJXdBW6witTCuL26DIlaz5C5oC/Ta8urVuL64DYpYzZK74Mah0Ezc\nBkWsZsldwMl71LyyrHOjRrpi7QrOseaZV5adqoQPHYF2nLzPNa8sO9kJFzoCFeaz9pq4DYpY\nzZIM/9zm72nFd+HkPWpeWTZXT5hYT+kwbX1xG7S25Y8ziNVLmFj74rH70ZyiOF+G/+75TNwG\nvRTrwi/Eagg8FHZP3ov64bvnmrgNiljNkgyINT77Rx+IdUmYWNcU7WfEWkA21+ZdiXU+h/pe\nrB8VvtON4c8M3IO05p4OhbVVy9xjnfmTPVbNPYlVP9ytWI0miNUsyTBMrGk/fOkV68z+8eum\nNojlPK88m6t7oFg2ly7oX8fh63ZJiOU8rzybqztarHJBh8KMUN+I1evXfcoyJpure6RYqTwd\nEevDYrp8TmXucNn/J79a90u69yvviNVd4ihW+nyol7af001/rv/g89x9lxsvw/wbdIhY9TNi\njTl5L78cuhCq7Pv6in6xCqtXiBWfzdUz+lBYtvdQl8fCUWJNZv4NiljdJRkmidVeiliI1cZD\nrN5D4xXdPzjsHlJ62Jk/lTX/BkWs7pIMk8Wqn0cfCvf1v6UorDdMmn+DIlZ3SYYpYrUvN4ze\nYz2n6oN++0fzXXTn36CI1V2SYZhYHtz4zPuC/pUOYnWXZECsAdnfZxCruyRDmFhLOhSOF6sB\nsWYXa0kn74jVn83VEybWki43IFZ/NldPnFhTmXGDIlZ/NlfP3xl8RUAsxJpFrKfzZa8HzrHc\n550rO00IL67vNnNayrtC/3nnyubq+d9ttHus4vybNj/WeR0LscLEWvcFUsQKE+up+s3m5WFn\n/mUFM25QxOrP5uoJE+vzAumHcX0zblC7WJd+zTXvXNlcPWFiNRdIzb9mbMYNilj92Vw9cWJN\nZY4NevLp39+I1Z/N1YNYORArm83Vg1g5ECubzdWDWDkQK5vN1TNMrNR5toBYiHUtVmo/2UAs\nxEKsGsTKZnP1DD0UpubRfFsQb7Hm4N8LcreIzLPBO0aOF8t6WxD2WOyx8nus8lqosu/rDlsW\nq0Y+78zZXD1Dxfq0ynxbEMRCrOFita1BLMT6wkesIbcFQSzE6hPreqc08rYgiIVYg8Uac1sQ\nxOq8ObxPWcZkc/UME8sDxEIsxKpBrGw2Vw9i9fK7BWLdyubqQaxeEGtINlcPYvWCWEOyuXoQ\nqxfEGpLN1cNNQXpBrCFZ50aNLFisAUYNEqtGMG9M1rlRI4iFWBIQC7EkIBZiSUAsxJKAWIgl\nAbEQSwJiIZYExEIsCYiFWBIQC7EkLEOs39cg1q2sc6NGEKuhXu99yjIm69yoEcRCLAmIhVgS\nEAuxJCAWYkn4XqziyJDnM5oNilgjsgJLDHwrVlE/fPdco9mgiDUiK7DEAGIhloRh51iINW7e\n0Ky3IjYcxfpRIRmyc9fRyTcf7aVeq+Yb2B4DxDqfnG9lj9X8s7JJ84Zm/SWxwKGwA2L5gFiI\nJYF3hYglAbEQSwJX3hFLAj8rRCwJiIVYEu5brB6hEOu7rHOjRhDrhlhD/EKs2yAWYklALMSS\ngFiIJQGxEEsCYiGWBMRCLAmIhVgSlijWaFkQa34QC7EkIFZerKxfiHUbxEIsCYiFWBIQC7Ek\nIBZiSbhTsfqNQqwhWedGjSDWjSxiTQOxEEsCYg0Sq98vxLoNYiGWBMRCLAmIhVgSEAuxJCAW\nYklALMSS4C2WE733hhTdIzJ758heorfOEmCPNWaP1dlxsce6DWIhlgTEQiwJ9yZW1ijEGpJ1\nbtQIYiGWBMRCLAmIlc9m/UKs2yAWYklALMSSsCixnGRBrBlALMSSgFh2sZrfFGYUQJV1btQI\nYiGWBMRCLAmIhVgS7kasVlFLE+v3EL8QaxrmjYRYXlnnRo0g1qAsYo0FsRBLAmK5iZX3C7Gm\nYd5IiOWVdW7UCGIhlgTE8hbrhl+INY3xG2lIURJZEEsJYo3JItZgEEskVtcvxJrG+I2EWOO3\nWTbr3KgRxNKK9SkYYk1j/Ebahlh11irLmKxzo0bixLq58e9YrBrE+p7pYhVHvr4avpGGiDWj\nLGOydrFqEOt7is+HE5sQq8Yu1iWI1cd4sUZs/PlliRDryi/EKlVihckyITtVrP4sYpU/Kqau\nDlZC3DnWoL99ZMdnpzbqA2KtLju1UR8Qa3XZqY36gFiry05t1AfEWl12aqM+xF15d9+gZOsl\nd0H8D6HJOmedGzWCWKvLOjdqBLFWl3Vu1AhirS7r3KgRxFpd1rlRI4i1uqxzo0YQa3VZ50aN\nyH8RpurzDqxXu96pIBbrlYBYrFcCYrFeCXf6y8Zh6SAWSEAskIBYIAGxQAJigQSJWO0PlTav\nLz9o6rveqSvurrco3ect2uv1nre77B5QiNX+GHzzuvPReLf1Tl1nd71X6/dab/PaeV7P7esK\nYs0nllf5mxXrhOob7+wB3LZlcfEkmdfp78HFehGrtcxjvc1rj1OW9no/T4W6/6+J6/1cpfe8\n9cNmxOr+zXc9DLRee69XOG9x+TBprYVqXjeWLtbVQof1KsW6/p9NW/HWxLrs3O8b7/XKeb2i\neftfTl3zxsTqnli4HQKuXnuvdwnzbvddYc+Jhcs33rO/cnmn1fPOzX3ellguq92iWEX7rZXj\nleH2ei+uZC9g3i8THK/oO8/rCj8rBAmIBRIQCyQgFkhALJCAWCABsUACYoEExAIJiHVBYns4\nwYa8ALG8YENegFhebHdDHtLD6fkhfZTvTykVu/Is1lmu0+PhOaXnQ+CQy2W7YpVPaX983B/9\neksndldiFdXih9gxF8qGxXqrTCp36e240/pVlh+NVF9ivVSJXXoNHnSRbFis8uH8eabq5f7t\n5fFarIfzq6fAGRfLlsV6Te/le3o5vno8Hwu7YqXULIexbHmjHdLz8UB3PDd/Tg+vb3vE8mTT\nG+057U/HufM7wAux9l+HQrCw6U33ftwbvZeVSu/l4fMcqzieyp+/2lUn77/SY/ScS2TTYh13\nSadrCbvUPsc6ffVSvTqcLjekj+gxl8i2xXqtrjOU1TExPb5/HgZ3xfGM/nxAPP1B6IhLZdti\ngQzEAgmIBRIQCyQgFkhALJCAWCABsUACYoEExAIJiAUSEAsk/Af0ls+gU/JpKgAAAABJRU5E\nrkJggg=="
	},
	"metadata": {}
	}
	]
	},
	{
	"metadata": {
	"ExecuteTime": {
	"start_time": "2019-04-18T04:28:46.526Z"
	},
	"trusted": true
	},
	"cell_type": "code",
	"source": "n <- 10^6\nN <- 1000\np <- 0.3\nt <- 2.0\n\nstart_time <- proc.time()\nMmean <- rep(0, n)\nFmean <- rep(0, n)\nfor (i in seq_len(n)) {\n M <- dqrnorm(N)\n F <- dqrnorm(N)\n f <- F[F > t]\n Mmean[i] <- mean(M[M > t])\n Fmean[i] <- mean(f[1:(length(f) * p)])\n}\nproc.time() - start_time",
	"execution_count": 5,
	"outputs": [
	{
	"output_type": "display_data",
	"data": {
	"text/plain": " user system elapsed \n 37.12 0.09 37.25 "
	},
	"metadata": {}
	}
	]
	},
	{
	"metadata": {
	"ExecuteTime": {
	"start_time": "2019-04-18T04:28:46.527Z"
	},
	"trusted": true
	},
	"cell_type": "code",
	"source": "options(repr.plot.width=5, repr.plot.height=3)\ndata.frame(Mmean, Fmean) %>%\n tidyr::gather() %>%\n ggplot(aes(x = value, fill = key)) +\n geom_histogram(position = \"identity\", bins=50, alpha = 0.7) +\n xlim(2, 3)",
	"execution_count": 6,
	"outputs": [
	{
	"output_type": "stream",
	"text": "Warning message:\n\"Removed 538 rows containing non-finite values (stat_bin).\"Warning message:\n\"Removed 4 rows containing missing values (geom_bar).\"",
	"name": "stderr"
	},
	{
	"output_type": "display_data",
	"data": {
	"text/plain": "plot without title",
	"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAFoCAMAAAC8KnXeAAAAUVBMVEUAAAAzMzNGzNBIztJJ\ntLVLtbdM0tZNTU1oaGh8fHyMjIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enr6+vw8PDy\n8vL0mZP2m5X6n5n///9X9oKHAAAACXBIWXMAABJ0AAASdAHeZh94AAAOHUlEQVR4nO2dbWOi\nuhZGM3PPOdjWvno9vcP//6FXFCsKSRPybCLNWh/E0mcy282aEKlDXQtggCtdAPxMEAtMQCww\nAbHAhBixmtPjgZgtQBslVu9N//DdFqDje7GaFrEgmfhTYf8MsSAGA7F+dajrhJWRJFbTRs9Y\n/71lvMcL0axhhXZkgFj3GK1MrKhFfI+0VdVF6xKruTwilm20KrGGeiGWbbQmsZqmv7Qed+Vd\n2qrqopWINQdpq6qLIpYXaauqiyKWF2mrqosilhdpq6qLIpYXaauqiyKWF2mrqosilhdpqzKj\nv4/8U66A5ChieZG2KjOKWAVArAUKSI4ilhdpqzKjiFUAxFqggOQoYnmRtiozilgF+NFi/R6C\nWIuCWJYFzI0ilhdpq+ZGEasciGVZwNwoYnmRtmpuFLHKgViWBcyNIpYXaavmRhGrHIhlWcDc\nKGJ5kbZqbhSxyoFYlgXMjSKWF2mr5kYRqxwVidVvlyxgbhSxvEhbNTeKWOVALMsC5kYRy4u0\nVXOjiFUOxLIsYG4UsbxIWzU3iljl+NG/QOCfKUoXVQnMWJYFzI0yY3mRtmpuFLHKgViWBcyN\nIpYXaavmRhGrHIhlWcDcKGJ5kbZqbhSxyoFYlgXMjSKWF2mr5kYRqxyIZVnA3ChieZG2am4U\nscrxM8X6PQFiLQpimRSQGUUsL9JWJUcRqziIZVJAZhSxvEhblRxFrOIglkkBmVHE8iJtVXIU\nsYqDWCYFZEYRy4u0VclRxCoOYpkUkBlFLC/SViVHEas4iGVSQGYUsbxIW5UcRaziIJZJAZlR\nxPIibVVyFLGKU59YEX4hVj6IZVJAZhSxvEhblRxFrOIglkkBmVHE8iJtVXIUsYoTL1ZzIGZ7\nQtqq5ChiFSdarKZ/+G7bI21VchSxioNYJgVkRhELsUyiiDUl1q8OdZ1pTN4jkjtGLkmSWA0z\n1jLRmmas7l0fp8KFolWJ1YFYy0RrEovF+4JRxEIsk2hNYnHlfcFoVWKlIW1VchSxioNYJgVk\nRhHLi7RVyVHEKg5imRSQGUUsL9JWJUcRqziIZVJAZhSxvEhblRxFrOIglkkBmVHE8iJtVXIU\nsYqDWCYFZEYRy4u0VclRxCoOYpkUkBlFLC/SViVHEas4iGVSQGYUsbxIWxUfnTIKsUqAWNoC\nNFHE8iJtVXwUsRDLpKuIhVgmXUUsxDLpKmIhlklXEQuxTLqKWIhl0lXEQiyTriIWYpl0FbEQ\ny6SriIVYJl1FLMQy6Wq8WCG/ECsfK7EKEXGPSO4YuQjMWNoCNFFmLC/SVsVHEQuxTLqKWGGx\n3GJLH8TSFqCJIpYXaavio4iFWCZdRawIsbburW33T8497du923T7+40QxNIWoImairV128Oz\nxh3Y9Ja1r+5ZLABiaQvQRC3F2h4deu7s2rqX9sM9HL58dO9iARBLW4AmaijWaYZqN8cD7x47\npz4OT5rQn5oDYmkL0EQNxXKHWap/0tEepqzH9s09qQVALG0BmqihWM+N27VDsQ6T166fxpQg\nlrYATdRyjfXanf/6U+GRN7dt9BoglrYATdT0XeGDe+3fGr4eF+4HyU5bKYilLUATNRXrwzX7\ndn+83NCt27spq1NNDGJpC9BEbS+QPndL9d3TYZ56P+/dyQVALG0Bmuiyn254l192bxFLXYAm\nuqxYD/r3hIilLkATXVIsZ7F0Ryx1AZrokmI1x8sPchBLW4Amuuyp0IQbsc6f12lyf3YkbVV8\nFLHuUazGDcgcV9qq+Chi3aNYLwOvXjLHlbYqPopY9yhWK/zoqrRV8VHEulOxZEhbFR9FrHsV\na9uwxkIsATcCbVm8L1ar1bCZR07EjUBN9qq9R9qq+ChihcX6n59/NQf+DIt3bQGa6M8T69Ht\nNeNKWxUf7XX56wRi3VJMrF3z4PtoTtOcLsd/tz0hbVV89FqsK78Qqy16KvQt3pv+4bttj7RV\n8VHEQiyTrv41BWJ9UUwsP81wi1i20brEOq2hvhfrV4e2ymj+DsCtIu/zVNhbtcIZ68TfzFj3\nKVb/cKdiTa7aEeuWOLEUP3yZ/NO7h/FNbRBrwWhpsdIcmmR6jP34dkmItWD0x4o18aOdOxcr\nYBRiXUgUy7XHM2J/WnTXW9eGTpfT33mduF/SfV95nyvWpF+IdX7ivh76vcOt8/oz/sbX2n0b\nKi8Caasiooh12REgZfHeXhy6Eqqd+nrEtFhNrlf3KFa/RayUU2E7nKGuz4VJYsmQtioiiliX\nHQGyxBruRSzEGqIQa/LUOOL2G/vtxrnNNvtTWdJWRUQR67IjQLZY/Tb5VLjr/y9Fk3vDJGmr\nIqKIddkRIEes4eWG5BnryXUf9Ns9ZN9FV9qqiChiXXYEiBNLgecz76v7XzqIddkRALFSo4h1\n2RGgmFjrOxV+nkCsy44AxcRa3+I9WaxJv5apdZlhQ4elmFjru9yAWKMdAcqJpULaqlAUsUY7\nAvwbQCsAYiHWImI9ni57bVhjGdS6zLCZR07E+G4zx70/+F1h1WL9x4/tjNWcftPmx3quYyHW\naEeAYmKt7wIpYo12BCgm1mP3m83b/Tb7lxVIWxWKItZoR4BiYn1dIP3IHFfaqlAUsUY7AhQT\n63yBNPvXjElbFYrOFevar2VqXWbY0GEpJ5YKaasmOfr05xOxRjsCINb3IJYnGjosiPU9iOWJ\nhg4LYn0PYnmiocMSJ5a72c4BsRBrLJYbbuaBWIiFWFcglicaOiyxp0J3fpx9WxDEQqygWHNv\nC2Illj1/rgjdfDTixqSlX8xSzJix2nm3BWHGYsaaEOvLqtm3BUEsxIoXaxBgxgqL1WNca2r0\nzsWKuS0IYiHWlFjjSSnxtiCIhVjRYqXcFgSxEGsklgLEOnK9hkesfFYo1ucAxEIsWVcRKxwN\nHRbECoBY4WjosCBWAMQKR0OHBbECIFY4Gjos3BQkAGKFo0ZHNBHEGiKvdWYUsbxIW3XNtVgR\nuiBWARALsUxALMQyAbEQywTEQiwTEAuxTEAsxDIBsRDLBMRCLBMQC7FMQCzEMgGxhvSjIlY+\naxLrcwxiIVZ+VxErKmp0RBNBLMQyAbEQy4R4sZoDMdsT0lb1IFZUVGhHBtFiNf3Dd9seaat6\nECsqKrQjA8RCLBPS1liIZfCy1MOq1MjDQKxfHdIie/5MMfcekZP0Y1oUXxsJYp0W5zXMWP3/\n/cmpNS/KjIVY+S9LPaxKjTwQawLEyod3hYhlAmIhlglceUcsE/hZIWKZsAqxJoxCLG/U6Igm\ngliIZQJiIZYJiOUX6zPCL8TygFiIZQJiIZYJiIVYJiAWYpmAWIhlAmIhlgmIhVgmrFesVFsQ\na1EQC7FMQKxvxQr6hVgeEAuxTEAsxDIBsRDLhPsWa9IoxApHjY5oIojljyJWBoiFWCYgFmKZ\ngFixYk37hVgeEAuxTEAsxDIBsRDLBMRCLBMQC7FMsBJLw+TNIS3uERm6ceQ1pVuyFpixEmes\nm5mLGcsDYiGWCYiFWCYgFmKZcKdiBY8pYgWjRkc0EcRCLBMQ69tosArE8oBYiGXCCsWS2IJY\nxiAWYpmAWIhlAmLliWVzH1zE8pLZKsSaHzU6ookgFmKZgFiIZcK9iTU4WIg1L2p0RBNBLIlY\nnxF+IZaCua26Q7GCfiGWB8RCLBMQSylW0C/EUjC3VYiVmkWsKBArNYtYUfwIsab9QiwFya2K\nOFgGtiCWFYiVGEWsOBDLTqwbvxBLQXKrECulW6Go0RFNBLHMxTr7hVgKklu1FrHOJIvVR2fa\nkhA1OqKJlBfLdwQQC7GONAcuX0nFWsyWhOhcsXoQK5bm6+HIjxerZ65Y1yCWn3Sx4o/A0rYs\nLtbIL8T6wkSsQrYs8ZnAMIMsYp0ffnWohoWVUn6NFfFvsLoop8ILiCWMItYFxBJGEesCYgmj\niHUBsYRRxBow98q7tqs/I4pYXqStqi6KWF6kraouilhepK2qLopYXqStqi6KWF6kraouilhe\npK2qLopYXqStqi6KWPHYfN6BUa2GzQWxVj4qYjGqyaiIxagmo1YvFtQFYoEJiAUmIBaYgFhg\nAmKBCaZiDT9Uen5+/UFT5ah5w96O2rTiWpvhqEYtyBtViqVYw4/Bn5/ffDReNGreiLejjkbX\njHp+nn/4bRqrBbHGo9qJJTr6tYt1xOb138wCooY2VxuDWjUT1tWwiNVaiqVZYJxrPS+Gbv+m\nrFG/BhQthioX6/Zfv+b1Ty5bZKOa1dpcP+QM2tgUq+NniDXamT2qnVjjvypr3ErFuj7qqtc/\n6ZV0VJNap59mDlynWLfLC83rH3ulHfXua63+XeHE8kLV1uvngndbE+/dxLUOxFLUWrFYzfDt\nlewC8XDUq6vZd17rRYX8d4VfBbbCYrXws0IwAbHABMQCExALTEAsMAGxwATEAhMQC0xALDAB\nsUY4eiKAJo5ALAU0cQRiKai7iXu3OW437qN9f3Su2bYnsU5yHR/3T8497QsWuU7qFqt9dLvD\n4+7g15s7sh2J1XS7N2XLXCGVi/XWmdRu3dth0npt24+zVBexnrvE1r0ULnR1VC5Wuzl9pKl7\nunt7fhiLtTk9eyxY4yqpXawX996+u+fDs4fTufBWLOfO+yGF2hu2d0+HE91hbf7kNi9vO8RS\nUX3DntzueJ47vQO8Emt3ORVCKtW37f0wG723nUrv7f5rjdUclvKnr7bd4v3VPZSuc21UL9Zh\nSjpeS9i64Rrr+NVz92x/vNzgPkqXuTYQ66W7ztB250T38P51Gtw2hxX96YR4/EbREtcIYoEJ\niAUmIBaYgFhgAmKBCYgFJiAWmIBYYAJigQmIBSYgFpiAWGDC/wGwkyHZzfyOSgAAAABJRU5E\nrkJggg=="
	},
	"metadata": {}
	}
	]
	},
	{
	"metadata": {
	"trusted": true
	},
	"cell_type": "code",
	"source": "simpi <- function (n) {\n s <- 0\n for (i in seq_len(n)) {\n x <- runif(2)\n s <- s + ifelse(x[1]^2 + x[2]^2 < 1.0, 1, 0)\n }\n 4*s/n\n}\n\nstart_time <- proc.time()\nsimpi(10^7)\nproc.time() - start_time",
	"execution_count": 17,
	"outputs": [
	{
	"output_type": "display_data",
	"data": {
	"text/html": "3.1403876",
	"text/markdown": "3.1403876",
	"text/latex": "3.1403876",
	"text/plain": "[1] 3.140388"
	},
	"metadata": {}
	},
	{
	"output_type": "display_data",
	"data": {
	"text/plain": " user system elapsed \n 31.83 0.01 31.85 "
	},
	"metadata": {}
	}
	]
	},
	{
	"metadata": {
	"trusted": true
	},
	"cell_type": "code",
	"source": "simpi_dqr <- function (n) {\n s <- 0\n for (i in seq_len(n)) {\n x <- dqrunif(2)\n s <- s + ifelse(x[1]^2 + x[2]^2 < 1.0, 1, 0)\n }\n 4*s/n\n}\n\nstart_time <- proc.time()\nsimpi_dqr(10^7)\nproc.time() - start_time",
	"execution_count": 18,
	"outputs": [
	{
	"output_type": "display_data",
	"data": {
	"text/html": "3.1424212",
	"text/markdown": "3.1424212",
	"text/latex": "3.1424212",
	"text/plain": "[1] 3.142421"
	},
	"metadata": {}
	},
	{
	"output_type": "display_data",
	"data": {
	"text/plain": " user system elapsed \n 22.68 0.00 22.72 "
	},
	"metadata": {}
	}
	]
	},
	{
	"metadata": {
	"trusted": true
	},
	"cell_type": "code",
	"source": "simpi_dqr1000 <- function (n) {\n s <- 0\n N = 1000\n m = n/N\n for (i in seq_len(m)) {\n x <- dqrunif(N)\n y <- dqrunif(N)\n for (j in seq_len(N)) {\n s <- s + ifelse(x[j]^2 + y[j]^2 < 1.0, 1, 0)\n }\n }\n 4*s/n\n}\n\nstart_time <- proc.time()\nsimpi_dqr1000(10^7)\nproc.time() - start_time",
	"execution_count": 25,
	"outputs": [
	{
	"output_type": "display_data",
	"data": {
	"text/html": "3.1411412",
	"text/markdown": "3.1411412",
	"text/latex": "3.1411412",
	"text/plain": "[1] 3.141141"
	},
	"metadata": {}
	},
	{
	"output_type": "display_data",
	"data": {
	"text/plain": " user system elapsed \n 15.62 0.00 15.62 "
	},
	"metadata": {}
	}
	]
	},
	{
	"metadata": {
	"trusted": true
	},
	"cell_type": "code",
	"source": "simpi_dqr_sum <- function (n) {\n x <- dqrunif(n)\n y <- dqrunif(n)\n 4*sum(x^2 + y^2 < 1.0)/n\n}\n\nstart_time <- proc.time()\nsimpi_dqr_sum(10^7)\nproc.time() - start_time",
	"execution_count": 39,
	"outputs": [
	{
	"output_type": "display_data",
	"data": {
	"text/html": "3.1406348",
	"text/markdown": "3.1406348",
	"text/latex": "3.1406348",
	"text/plain": "[1] 3.140635"
	},
	"metadata": {}
	},
	{
	"output_type": "display_data",
	"data": {
	"text/plain": " user system elapsed \n 0.19 0.16 0.41 "
	},
	"metadata": {}
	}
	]
	},
	{
	"metadata": {
	"trusted": true
	},
	"cell_type": "code",
	"source": "start_time <- proc.time()\nsimpi_dqr_mean(10^8)\nproc.time() - start_time",
	"execution_count": 40,
	"outputs": [
	{
	"output_type": "display_data",
	"data": {
	"text/html": "3.14144876",
	"text/markdown": "3.14144876",
	"text/latex": "3.14144876",
	"text/plain": "[1] 3.141449"
	},
	"metadata": {}
	},
	{
	"output_type": "display_data",
	"data": {
	"text/plain": " user system elapsed \n 1.95 2.06 4.02 "
	},
	"metadata": {}
	}
	]
	},
	{
	"metadata": {
	"trusted": true
	},
	"cell_type": "code",
	"source": "start_time <- proc.time()\nsimpi_dqr_mean(10^9)\nproc.time() - start_time",
	"execution_count": 41,
	"outputs": [
	{
	"output_type": "error",
	"ename": "ERROR",
	"evalue": "Error: サイズ 7.5 Gb のベクトルを割り当てることができません \n",
	"traceback": [
	"Error: サイズ 7.5 Gb のベクトルを割り当てることができません \nTraceback:\n",
	"1. simpi_dqr_mean(10^9)",
	"2. mean(x^2 + y^2 < 1) # at line 4 of file <text>"
	]
	}
	]
	},
	{
	"metadata": {
	"trusted": true
	},
	"cell_type": "code",
	"source": "",
	"execution_count": null,
	"outputs": []
	}
	],
	"metadata": {
	"kernelspec": {
	"name": "ir",
	"display_name": "R",
	"language": "R"
	},
	"language_info": {
	"name": "R",
	"codemirror_mode": "r",
	"pygments_lexer": "r",
	"mimetype": "text/x-r-source",
	"file_extension": ".r",
	"version": "3.5.3"
	},
	"toc": {
	"nav_menu": {},
	"number_sections": true,
	"sideBar": true,
	"skip_h1_title": false,
	"title_cell": "Table of Contents",
	"title_sidebar": "Contents",
	"toc_cell": false,
	"toc_position": {},
	"toc_section_display": true,
	"toc_window_display": false
	},
	"gist": {
	"id": "9b2f093157d8ee9e4587fbf06eb0ba5e",
	"data": {
	"description": "MFmean by R.ipynb",
	"public": true
	}
	},
	"_draft": {
	"nbviewer_url": "https://gist.github.com/9b2f093157d8ee9e4587fbf06eb0ba5e"
	}
	},
	"nbformat": 4,
	"nbformat_minor": 2
	}